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ABSTRACT 


This  report  describes  computer  programs  developed  for  the  analysis  of  heated 
beams,  plates,  and  stiffened  cylindrical  shells.  The  matrix  displacement  approach 
to  structural  analysis,  which  forms  the  theoretical  basis  of  these  programs,  is 
developed  in  detail.  Derivation  of  new  relationships  employed  in  these  programs  is 
also  detailed.  The  capabilities  and  limitations  of  the  respective  programs  are  out¬ 
lined  and  illustrative  applications  are  presented. 
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CHAPTER  I 
INTRODUCTION 


The  "Study  of  Thermal  Stress  Determination  Techniques  for  Supersonic 
Transport  Aircraft  Structures"  has  consisted  of  three  relatively  independent  efforts, 
directed  towards  the  development  of 

(1)  An  annotated  bibliography  of  literature  pertinent  to  thermal  stress 
analysis  and  related  topics. 

(2)  Parametrically-presented  design  data  for  heated  sandwich  panels  and 
cylinders. 

(3)  Three  POUT  RAN -coded  computer  programs  for  the  solution  of  heated 
beam,  plate,  and  cylindrical  shell  problems. 

Items  (1)  and  (2)  a  re  presented  in  references  (1)  and  (2),  respectively.  A  portion  of 
item  (3),  in  the  form  of  a  verbal  description  of  the  coded  computer  programs,  is 
presented  in  this  report. 

Each  of  the  three  programs  described  in  the  present  report  is  based  on  the 
"matrix  displacement"  or  "stiffness"  method  for  the  analysis  of  structures  which 
arc  idealized  as  systems  of  connected  discrete  elements.  Programs  described  in 
the  present  report  are  available  to  participants  in  structural  design  activities 
related  to  the  Supersonic  Transport  (SST)  and  to  all  others  who  are  designated  as 
being  eligible  by  the  Flight  Dynamics  Laboratory,  Aeronautical  Systems  Division, 
USAF.  The  programs  will  be  transmitted  to  eligible  recipients  by  the  latter 
agency.  A  transmitted  program  consists  of  punched  cards,  listings,  detailed 
instructions  with  respect  to  input  and  output,  and  other  information  needed  to  make 
the  program  operative  at  a  facility  that  will  accept  a  program  coded  in  conformity 
with  the  FORTRAN  II  Monitor  System. 

Many  references  have  detailed  the  basis  of  the  displacement  method  as  it 
applies  to  the  linear  analysis  of  u cheated  elastic  structures.  (See,  for  example, 
references  3,  4,  or  5).  Phenomena  which  are  not  often  considered  in  routine 
structural  analysis  are  treated  by  the  subject  programs,  however,  and  the  approach 
to  the  analysis  of  these  special  phenomena  by  the  matrix  displacement  method  has 
not  been  described  in  any  single  reference.  Hence,  in  the  next  chapter,  the  method 
is  developed  from  fundamental  principles  and  to  an  extent  that  includes  all  pertinent 
special  phenomena,  such  as  instability,  thermal  stress,  and  inelastic  behavior. 

The  accuracy  and  efficiency  of  any  solution  performed  by  use  of  any  computer 
program  for  matrix  structural  analysis  is  largely  dependent  upon  the  suitability  of 
discrete  element  force-displacement  equations  employed.  The  techniques  used  in 
derivation  of  the  element  force-displacement  equations  contained  in  the  subject 

Manuscript  released  by  authors  in  Jan.  1964  for  publication  as  an  ASD 
Technical  Documentary  Report. 
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programs  are  not  well  known  and  some  of  these  techniques  have,  in  fact,  been  for¬ 
mulated  specifically  for  the  requirements  of  this  study.  Consequently,  in  Chapter 
III  a  complete  and  detailed  development  of  the  techniques  is  presented.  Chapter  in 
also  demonstrates  how  the  techniques  were  applied  to  the  derivation  of  element 
relationships  appearing  in  the  respective  computer  programs. 

The  three  coded  programs  are  as  follows: 

1.  A  program  for  the  analysis  of  one-dimensional  structural  components. 

2.  A  program  for  the  analysis  of  flat  plates, 

3.  A  program  for  cylinder  analysis. 

Program  (1)  is  described  in  Chapter  IV.  The  objective  of  this  program  is  to 
permit  analyses  of  beam-type  structures,  i.e.,  structures  whose  cross-sectional 
dimensions  are  small  with  respect  to  their  length  and  whose  behavior  is  governed 
by  the  elementary  concepts  of  beam  flexure. 

It  often  proves  feasible  in  the  analysis  of  airframes  to  isolate  portions  of 
major  components  and  treat  these  as  one-dimensional  elements;  i.e.,  as  beams  or 
beam-columns.  The  most  common  examples  in  modern  constructional  forms  are 
spar  caps,  stringers,  and  longerons,  but  such  idealizations  may  also  be  admissable 
in  connection  with  trussed  internal  members  (spars,  ribs)  and  control  surfaces. 

The  objectives  of  the  program  described  herein  pertain  to  these  types  of  elements. 

With  this  program  it  is  possible  to  analyze  beams  of  nonuniform  section  over 
many  supports  for  thermal  stress,  inelastic  behavior,  instability  and  other  types  of 
structural  behavior.  Chapter  IV  provides  a  detailed  picture  of  these  capabilities 
and  also  presents  illustrative  examples.  Certain  of  the  examples  involve  simple 
conditions  with  known  solutions;  these  are  performed  to  demonstrate  the  accuracy 
of  the  program.  Another  example,  for  which  there  is  no  complete  alternative  solu¬ 
tion,  is  performed  to  demonstrate  the  capability  of  the  program  to  deal  with  com¬ 
plex  conditions. 

Descriptions  of  Programs  (2)  and  (3),  which  are  presented  in  Chapters  V  and 
VI,  are  patterned  after  the  description  of  Program  (1).  Program  (2)  is  capable  of 
performing  analyses  of  irregularly  shaped  stiffened  plates  of  nonuniform  thicknesses 
for  stresses  and  displacements  due  to  applied  loads,  temperature  gradients,  and 
time-independent  inelastic  behavior,  and  for  the  prediction  of  buckling  stresses. 

The  use  of  the  discrete  element  approach,  rather  than  design  charts,  may  be 
necessary  even  for  isotropic  rectangular  plates  of  constant  thickness  where  the 
temperature  profile  (and  therefore  the  resulting  stress  distribution)  is  of  a  highly 
irregular  form.  Also,  if  the  loadings  on  the  plate  have  been  developed  through  a 
matrix  structural  analysis,  wherein  the  plate  was  employed  as  a  single  discrete 
element  in  a  major  component  of  the  airframe  (e.g.,  a  wing  or  fuselage),  the  edge 
loadings  will  in  general  be  nonuniform  and  it  again  may  be  desirable  to  utilize  the 
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discrete  element  approach  in  the  prediction  of  instability.  The  program  is,  of 
course,  not  limited  in  applicability  to  skin  panels;  it  is  directly  useful  for  analyses 
of  all  types  of  planar  structures  —  stiffened  bulkheads,  plane  trusses,  beam  grid- 
works,  etc. 

Program  (3),  the  cylinder  analysis  program,  can  be  employed  to  predict  the 
stresses  and  displacements  for  heated  cylindrical  shells.  The  structure  analyzed 
with  this  program  need  not  be  perfectly  cylindrical.  They  may  possess  orthotropic 
skins  and  can  be  ring  and  longitudinally-stiffened.  The  temperature  dependence  of 
material  properties  can  be  taken  into  account. 
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CHAPTER  II 

PERTINENT  CONCEPTS  OF  THE  MATRIX 
DISPLACEMENT  METHOD 


A.  GENERAL  THEORY  FOR  LINEAR  ELASTIC  UNHEATED  SYSTEMS 


Most  treatments  of  the  displacement  approach  to  the  analysis  of  discrete 
element  systems  make  extensive  use  of  matrix  notation  and  the  concepts  of  matrix 
algebra.  Matrices  present  an  especially  concise  and  convenient  means  of 
expressing  algebraic  procedures  and  are,  in  addition,  the  natural  mathematical 
language  of  the  digital  computer.  Because  of  this  widespread  use  of  matrix  notation 
in  the  formulation  of  structural  analysis  programs,  the  topic  of  present  interest  is 
often  referred  to  as  "matrix  structural  analysis".  For  simplicity,  the  initial  por¬ 
tion  of  the  following  development  of  the  displacement  method  will  be  presented  in 
detailed  algebraic  form;  then,  the  formulation  will  be  summarized  with  use  of 
matrix  algebra. 


As  already  noted,  this  report  is  concerned  with  applications  based  on  discrete 
element  idealizations  of  the  structures  to  be  analyzed.  Such  discrete  elements  are 
usually  defined  by  boundary  or  corner  points  which  are  sufficient  in  number  to 
characterize  the  stress  and  deformational  behavior  of  the  element.  A  hypotheti¬ 
cal  discrete  element,  with  four  boundary  points,  is  sketched  below.  For  convenience, 
this  section  deals  with  relationships  in  two  dimensions;  in  all  cases,  however,  the 
approach  to  three-dimensional  problems  is  readily  apparent. 


Relationships  between  the  forces  acting  at  the  corner  points  and  the  dis¬ 
placements  of  the  corner  points  are  derived  in  detail  in  Chapter  III.  Generally, 
these  relationships  are  first  derived  with  reference  to  axes  which  are  most  con¬ 
venient  to  the  element  itself.  Then,  by  use  of  direction  cosines,  the  relationsips 
are  expressed  as  if  the  element  were  arbitrarily  oriented  in  the  complete  structure. 
When  this  transformation  has  been  accomplished,  the  equation  for  one  of  the  forces 
shown  in  the  sketch  below,  F  for  example,  would  have  the  form 


F\3  k(x3)(ul)  U1  '  k(x3)(«2>  2  *  .  k(x3)(v4)V4 


(II-l) 
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The  coefficients  ai'e  functions  of  the  geometry  and  material  properties  of 

the  element  and  are  known  as  element  stiffness  coefficients.  In  terms  of  physical 
meaning,  a  stiffness  coefficient  kjj  is  the  force  at  i  necessary  to  produce  a  unit 
displacement  of  the  degree  of  freedom  j.  The  corner  point  forces  may  either  be 
actual  concentrated  forces  or  the  static  equivalent  of  stresses  acting  upon  the  area 
subtended  by  the  point;  normally,  they  are  the  latter. 


A  complete  set  of  force-displacement  relationships  for  an  element  with  n 
node  points  in  two  dimensions  will  appear  as: 


F(.xl)  k(xl)(ul)  V  • 


k(xl)  (un)  Un  4  k(xl)(vl)  2  + 


k(xl)(vn)Vn 


F  ,  =  k,  , vi.  +  ....  k 

(xn)  (xn)  (ul)  1 

F,  =  k,  u  +  .  .  .  .  k 

(yl)  (yl  (ul)  1 

F  =  k  u  +  .  .  .  .  k 

(yn)  (yn)  (ul)  1 


+  k 
+  k 
+  k 


+  • 


k 

k 

k 


(II- 2) 


All  degrees  of  freedom  at  the  boundary  points  appear  in  Equation  (II— 2) ,  i.e.,  the 
element  is  not  fixed  against  displacement  as  a  rigid  body. 


Once  the  element  force-displacement  relationships  (Equation  II— 2)  have  been 
numerically  evaluated  for  each  element  of  the  structure,  they  can  be  algebraically 
combined  in  a  manner  dictated  by  the  requirements  of  juncture  point  equilibrium 
and  compatibility'.  These  operations  produce  a  set  of  force-displacement  equations 
for  the  elemen*  juncture  points  of  the  assembled  structure.  To  illustrate  how  this 
is  accomplished,  consider  the  development  of  the  force-displacement  equation  at 
point  i  in  the  x-direction  in  the  assembled  analytical  model  sketched  below. 
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The  theoretical  basis  for  the  establishment  of  the  desired  relationships  is 
the  explicit  requirement  of  juncture  point  equilibrium.  This  condition  states  that  the 
applied  load  (P^)  is  equal  to  the  sum  of  the  internal  forces  acting  upon  the  respec¬ 
tive  elements  common  to  the  point 


=  F. 


Fxi 


+  F 


xx 


(II- 3) 


A 

whei'e  F  .  is  the  x-direction  (internal)  force  on  element  A,  etc.  The  force-dis¬ 
placement  equations  for  each  of  the  elements  will  have  been  previously  evaluated  so 
that  expressions  for  F  etc.  in  terms  of  the  displacements  are  available. 

Substitution  of  these  into  Equation  11-3  yields: 


P  .  =  (k,  ..  u.  i . 

XI  (XI)  (UJ)  j 

(xi)(ui,Ui  + 

■  • ) 

n 

+  fls.  Ur  t . 

,..kB 

(xi)  (uf)  f 

(XX)  (ui)  u .  .  . 

. .) 

.  C 

+  (k  .  .  ,  ii  t . 

(XI)  (ugj  g 

■  ■  ■  kC(xi)(ui)  Ui  +  - 

. .) 

(II-3a) 

and,  since  the  displacement  u.  is  the  same  for  A,  B,  and  C  at  point  i  (the  condition 
of  compatibility)  we  have: 


P  .  - 
xi 


,  A  t 

k  ......  «j  t  ....  (k 

(xi)  (uj) 


(xi)  (ui) 


+  k 


B 


(xi)  (ui) 


+  k 


)  u 

(xi)  (ui)'  i 


(II-3b) 


This  is  the  final  form  of  the  desired  equation. 


It  is  important  to  note  that  each  of  the  three  elements  meeting  at  the  indicated 

juncture  point  possess  stiffness  coefficients  with  common  subscripts  (e.g.k^ 

B  b  (xi)(ui)’ 

k  ,,,.).  These  correspond  to  the  common  juncture  of  certain  points  on  the 
(•xi)  (ui)'  H  ^ 

element.  Point  i  is  of  course  a  common  juncture  of  all  three  elements,  but  two 
elements  also  meet  at  points  h  and  1  .  Some  coefficients  associated  with  an  ele¬ 
ment  will  not  have  a  counterpart  coefficient  in  the  relationships  for  the  other  ele¬ 
ments.  These  will  pei'tain  to  points  such  as  g,  f,  and  j  which  are  connected  to  only 
one  of  the  elements  meeting  at  i  . 


In  view  of  the  above  reasoning,  the  following  "automatic"  approach  to  calcula¬ 
ting  the  applied  load-versus-displacement  equations  for  the  complete  structure 
suggests  itself: 


(1)  Each  element  stiffness  coefficient  is  assigned  a  double  subsci'ipt — the 
first  is  the  force  to  which  it  is  equated  and  the  second  is  the  displace¬ 
ment  it  multiplies. 
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(2)  Provision  is  made  for  an  equation  for  each  force  in  every  degree  of 
freedom  in  the  complete  system,  and  for  the  possibility  that  each  force 
will  be  related  to  every  displacement  in  the  system.  The  result  is  a 
rectangular  array  of  spaces,  each  designated  by  two  subscripts — the  first 
pertains  to  the  force  equation,  the  second  to  the  displacement  in  question. 

(3)  The  numerical  formulation  of  the  equations  can  begin  with  point  1  in  the 
x-direction.  First,  a  search  is  made  through  the  list  of  elements,  which 
are  designated  by  their  corner  points.  When  an  element  is  reached  whose 
designation  contains  a  1,  the  F  ^  equation  is  selected  and  each  coefficient 
of  the  equation  is  placed  in  the  space  reserved  for  its  second  subscript 
(each  of  die  first  subscripts  is,  of  course,  xl). 

(4)  The  procedure  of  item  (3)  is  continued  for  point  1  in  the  x-direction  until 
the  list  of  elements  is  exhausted.  Each  time  the  stiffness  coefficient  from 
a  subsequent  element  is  placed  in  a  location  where  a  value  has  already 
been  placed,  it  is  added  to  that  value.  The  only  locations  that  will  be 
occupied  as  a  result  of  these  operations  are  those  relating  to  the  point  in 
question  and  points  that  exist  on  elements  meeting  at  the  point  in  ques¬ 
tion.  Thus,  if  the  structure  is  large  with  many  clement  juncture  points, 
there  will  be  many  zeros  in  each  equation.  This  is  an  advantageous 
feature  in  terms  of  the  effort  needed  to  solve  the  complete  set  of  equa¬ 
tions. 

(5)  The  process  of  steps  (3)  and  (4)  is  repeated  for  all  other  points  in  the 
x-direction  and  then  for  each  of  the  other  directions.  The  result  will  be  a 
complete  set  of  equations  for  the  entire  structure,  but  with  no  recognition 
of  support  conditions. 

(G)  The  support  conditions  are  accounted  for  by  first  noting  which  displace¬ 
ments  are  zero  and  then  removing  the  stiffness  coefficients  multiplying 
these  displacements  from  the  equations,  resulting  in  more  equations  than 
unknowns.  To  provide  for  an  equality  of  equations  and  unknowns  the 
equations  which  pertain  to  the  external  loads  (reactions)  at  the  support 
points,  are  removed.  The  general  solution  to  the  remaining  set  of  equa¬ 
tions  gives  the  displacement  influence  coefficients;  multiplication  of  the 
general  solution  by  specific  values  for  the  loadings  yields  specific  values 
for  the  displacements. 

(7)  By  substituting  the  displacement  values  back  into  the  element  force-dis¬ 
placement  equations  the  internal  forces  acting  on  the  element- the  corner 
point  forces  -  can  be  determined.  These  may  require  retransformation 
from  "system"  to  "local"  coordinates,  and  finally  a  transformation  into 

stress  (  cr  ,  c  ,  r  ,  etc.), 
x  y  xy 
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B.  MATRIX  FORMULATION 


The  foregoing  can  now  be  reviewed  in  the  interests  of  obtaining  a  matrix  for¬ 
mulation  and  a  clearer  picture  of  the  computer  operations  in  a  practical  program. 
First,  any  complete  set  of  element  force-displacement  relationships  (Equation  II— 2) 
can  be  written  in  matrix  notation  as 


W  M  kl 


(II-2a) 


where  [kj  ,  the  "element  stiffness  matrix",  is  composed  entirely  of  the  element 


stiffness  coefficients  k 


,  etc.  The  subscript 


- - — —  •'(xl)(ul)’ . 

quantities  refer  to  the  node  points  of  the  element. 


denotes  that  the  indicated 


Once  the  element  relationships  have  been  evaluated,  the  elements  are 
assembled  to  form  the  complete  analytical  model  of  the  structure  by  joining  all  ele¬ 
ments  at  their  respective  juncture  points  and  applying  in  the  process  the  require¬ 
ments  of  juncture  [xiint  equilibrium  and  compatibility.  Thus,  the  components  of 
internal  loads  j  F  j  and  external  loads  j  P  j  at  each  point  arc  related  by  equilib¬ 
rium  requirements;  i.e.,  7  F^  =  P  ,  etc.  The  respective  coordinate  displace¬ 

ments  of  the  corner  points  of  all  elements  meeting  at  a  point  are  equal,  a  require¬ 
ment  that  satisfies  compatibility.  It  follows  that  the  stiffness  matrix  [kJ  for  the 
complete  structure  can  lie  assembled  by  merely  adding  element  stiffness  coeffi¬ 
cients  having  identical  subscripts.  This  results  in  a  set  of  equations: 


M  W 


(R-4a) 


The  matrix  J"k1  will  henceforth  be  referred  to  as  the  "master"  stiffness  matrix. 
Displacement  boundary  conditions  can  be  readily  imposed  by  assigning  the  pertinent 
A's  their  known  values  (usually  zero).  The  matrix  £l<j  will  be  altered  in  the 
process,  and,  taking  note  of  this  by  utilizing  the  subscript  R,  the  solution  to  the 
altered  Equation  (II-4a)  becomes  (if  matrix  inversion  is  utilized). 

K)  [K,J  {pR}  H  {pR}  <n-5» 


where 


represents  the  set  of  displacement  influence  coefficients. 


The  subject  computer  programs,  in  their  present  form,  are  restricted  to  the 
use  of  a  matrix  inversion  procedure  for  the  solution  of  Equation  (II— 5) .  One 
advantage  of  the  use  of  matrix  inversion  is  that  the  analysis  for  load  conditions 
additional  to  the  first  is  accomplished  at  small  additional  expense.  Thus,  the 
respective  programs  each  allow  the  solution  for  many  load  conditions  in  one  com¬ 
putational  cycle  when  linear  elastic  analyses  for  applied  load  are  attempted.  It  is 
to  be  recognized,  however,  that  the  direct  solution  of  Equation  II— 5  (e.g.,  through 
Gaussian  elimination  procedures  or  by  iterative  techniques)  may  prove  more 
advantageous  under  many  circumstances.  In  such  cases  it  is  possible  to  replace 
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the  inversion  subroutine  used  in  the  program  with  whatever  direct  solution  sub¬ 
routine  may  be  available. 


To  obtain  the  stresses  from  the  displacement  solution  there  can  first  be 
selected,  from  the  total  column  of  displacements,  the  displacement  vectors  for  the 
respective  elements  })  *  Then,  each  vector  can  be  multiplied  into  the 

stiffness  matrix  (Equation  ll-2a)  to  determine  the  node  point  forces  and 

in  an  additional  step,  the  node  point  forces  can  be  transformed  into  the  correspon¬ 
ding  stresses.  It  is  believed  to  be  more  efficient,  however,  to  form,  at  the  outset, 
direct  relationships  between  the  element  stresses  and  the  node  point  displacements, 
as  follows 


W  -  M  k) 


(II- 6) 


where 


k} 


are  the  stress  values  which  describe  the  distribution  of  stress 

within  an  element  and  ]  is  known  as  the  "element  stress  matrix".  The  proce¬ 
dure  followed,  therefore,  is  to  establish  first  the  stress  matrices  at  the  start  of  a 


computation.  When  the  displacement  vectors  for  the  respective  elements 


(kl) 


are  evaluated,  they  are  premultiplied  by  the  corresponding  clement  stress  matrices 
to  obtain  the  solutions  for  stress. 
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C.  INITIAL  STRAIN  PROBLEMS 

Problems  associated  with  thermal  and  plastic  strain  can  be  classed  as  "initial 
strain"  problems  since  they  can  exist  prior  to  the  imposition  (or,  in  the  case  of 
plasticity,  reimposition)  of  applied  load.  If  a  discrete  element  is  in  the  state  of 
initial  strain,  the  element  force-displacement  relationships  take  the  following  form 


k}  Mk}-k} 


(II- 7) 


where  {  Fe*  }  signifies  the  system  of  "initial  forces"  at  the  element  node  points. 
Physically,  these  values  represent  the  forces  required  to  impose,  at  the  node  points, 
displacements  which  are  equal  and  opposite  to  those  accruing  from  the  initial  strains. 
In  other  words,  they  are  the  forces  required  to  suppress  the  node  point  displace¬ 
ments  due  to  initial  strain.  A  procedure  for  the  derivation  of  element  initial  forces 
will  be  presented  in  the  next  chapter. 

Upon  assembly  of  Equations  (II- 7)  to  form  the  master  set  of  stiffness  equations, 
and  reduction  in  cognizance  of  boundary  conditions,  there  is  obtained 

{  P}  [k]  {a}  -  {P1}  (II- 8) 

where  now  the  values  •[  P‘|  arc  the  "net  thermal  forces  at  the  node  points.  The 
solution  to  (II-8)  is  given  by 


{A}  [k]-‘{{pMp‘}} 


(H-9) 


To  obtain  the  solution  for  stress  the  adopted  approach  is  to  formulate  stress- 
displacement  equations,  i.e.  "element  stress  matrices".  In  the  presence  of  initial 
strain  these  take  the  form 


k}-M  k)  -  K' 


(III- 10) 


where  {%'}  represents  the  stresses  required  to  obviate  the  initial  strains. 
D.  ELASTIC  INSTABILITY  ANALYSIS 


T 


I 


I 

r 


The  concepts  of  clastic  instability  pertain  to  conditions  in  prismatic  or  thin 
walled  structures,  where  the  behavior  across  the  thickness  can  be  subdivided  into 
"flexural"  and  "midplane"  behavior.  By  virtue  of  displacements  normal  to  the  mid¬ 
plane  the  midplanc  forces  have  components  which  tend  to  enhance  these  displace¬ 
ments.  When  their  magnitude  is  sufficiently  large  they  produce  infinitely  large  dis¬ 
placements,  whatever  the  magnitude  of  the  loads  applied  normal  to  the  midplane. 

The  values  of  midplane  load  which  cause  this  elastic  instability  are  the  "critical  loads." 

In  two  of  the  developed  computer  programs  -  the  one-ditnensional  and  plate 
programs  -  it  is  possible  to  separate  completely  the  midplane  and  out-of-plane 
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(flexural)  behaviors.  When  considering  elastic  instability,  the  determination  of  the 
midplane  stresses  and  displacements  is  unaffected  with  respect  to  the  procedures 
discussed  previously.  However,  element  stiffness  for  out-of-plane  behavior,  [kez]  ’ 
now  becomes  the  sum  of  two  component  stiffnesses: 

[  \  ]  '  [  S  ]  '  t "  ]  <n-U) 

where  [  kCf  ]  is  the  stiffness  for  conventional  flexural  behavior  and  [  n  ]  repre¬ 
sents  the  effects  of  the  midplane  forces  throughout  the  element  on  the  out-of-plane 
behavior.  The  terms  of  [  n  J  consist  of  the  dimensions  of  the  element  and  the  values 
of  midplane  force,  as  determined  in  the  (independent)  midplane  analysis.  Material 
properties  do  not  appear  in  the  [  n  J  matrix.  Techniques  for  formulating  Equation 
(II— 1 1)  for  beam  and  plate  elements  were  delineated  in  Reference  6  and  are  discus¬ 
sed  in  the  next  chapter. 


Due  to  the  segregation  of  midplane  and  out-of-plane  behavior,  two  separate 
sets  of  master  stiffness  equations  would  be  formed  in  an  instability  problem.  The 
midplane  equations  would  appear  as  (excluding  initial  strain  effects) 


j  p  l 

= 

r  k  i 

j  a  1 

l  xy  j 

L  *yj 

1  ^  xy  J 

(11-12) 


while  the  out-of-plane  equations  would  take  the  form 

K }  [[  kJ  *  M]  {a,}  di-13) 


To  solve  Equation  (11-13)  it  is  of  course  first  necessary  to  solve  Equation 
(11-12)  and  determine  the  associated  midplane  forces  so  that  the  matrix  [n^J  can  be 
formed.  In  the  form  indicated,  Equation  (11—13)  can  be  solved,  under  certain  condi¬ 
tions,  to  obtain  an  "equilibrium"  solution  for  out-of-plane  behavior  in  the  presence 
of  given  midplane  forces.  These  conditions  dictate  that  the  midplane  forces  be  of 
less  than  cx'itical  value.  When  they  are  of  critical  value  or  greater,  the  matrix  to  be 
inverted  will  be  singular. 


In  practice,  the  solution  of  these  equations  as  an  instability  problem  involves 
the  determination  of  the  value  of  the  midplane  forces  to  cause  instability.  It  is 
assumed  that  all  midplane  forces  arc  at  a  fixed  ratio  to  one  another  at  all  levels  of 
applied  load,  from  the  onset  of  loading  to  the  achievement  of  instability.  (Corres¬ 
pondingly,  the  shape  of  a  midplane  temperature  distribution  causing  midplane  forces 
remains  constant  up  through  instability).  Thus,  the  midplane  analysis  is  performed 
for  any  convenient  magnitude  of  the  applied  loads  and  it  is  assumed  that  at  instability 
the  actual  intensity  is  a  scalar,  X  ,  times  such  magnitude.  Equation  (11-13)  can  then 
be  written  as 

{P,.}-[K,.]  {  A„}*A[n]  {aJ  (11-14) 
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Also,  it  has  been  stated  that  elastic  instability  is  independent  of  the  value  of  the 
applied  midplane  loads.  Hence,  setting  j  \  =  0 


0  = 

K1 W 

*XIN]  W 

(11-15) 

T( 

>„}  -  Is 

(11-16) 

Using  matrix  iteration,  the  above  can  be  solved  for  the  eigenvalues  /Xj  and  the 
associated  eigenvectors  {A  z|.  ,  There  will  be  as  many  such  eigenvalues  as  there 
are  equations  in  (11—15) ,  but  the  only  eigenvalue  of  interest  is  the  largest  value  of 
representing  the  smallest  X^  and  therefore  the  lowest  magnitude  of  midplane 
load  at  which  elastic  instability  will  be  experienced. 
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CHAPTER  III 

DISCRETE  ELEMENT  FORCE-DISPLACEMENT 
RELATIONSHIPS 

A.  DERIVATION  PROCEDURES 
1.  Lineal'  Elastic  Stiffness 

The  force-displacement  properties  of  discrete  elements,  of  the  form 
required  for  use  in  matrix  displacement  analyses,  can  be  formulated  by  application 
of  one  of  three  general  approaches.  These  approaches  are  outlined  in  Reference  7 
and  developed  in  detail  in  Reference  5  .  In  the  derivation  of  relationships  for  the 
subject  group  of  three  matrix  displacement  computer  programs  it  proved  convenient 
and  sufficiently  accurate  to  employ  only  one  of  these  three  approaches.  The  selected 
approach  is  based  on  Castigliano’s  First  Theorem,  Part  I,  and  is  formulated  from 
fundamental  principles  in  this  chapter  of  the  report. 

To  formulate  this  approach  for  the  case  of  linear  elastic  stiffness  alone 
(initial  strain  and  instability  effects  being  temporarily  disregarded — those  are 
examined  in  later  sections)  consider  the  discrete  element  free  body  diagram, 
sketched  below  (Figure  III—  1) .  The  element  is 


F, 
i 

Figure  III—  1 .  Discrete  Element 

subjected  to  the  indicated  node  point  loads  F^,  F£,  .  .  .  Fj,  .  .  .  which  include 
both  applied  and  leactive  forces.  The  corresponding  node  point  displacements  are 
Af,  Ag,  ...  A.  ,  .  .  An  amount  of  elastic  strain  energy  (U)  is  stored  in 

the  structure  as  a  consequence  of  the  loading  and  displacements. 

If  the  structure  is  restrained  against  displacements  at  all  points  of 
load  application  except  at  and  in  the  direction  of  the  i  th  load,  Fj,  an  infinitesmal 
increase  in  the  load  FjtoFj+  will  result  in  an  incremental  displacement 
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SAj.  In  accordance  with  the  principle  of  conservation  of  energy,  the  change  in 
external  work  (  SW)  done  by  the  infinites mal  force  change  must  equal  the  change 
in  the  stored  elastic  strain  energy,  S  U.  Hence, 

SW  s  u 

and,  since  SW  (Fj  +  SFi)  SAj  ~^Fj  )  ^  SA  ^.Equation  III 
(Fj  )  (sAj)  SW  SU 
when  S  A  — 0 

dW  d  U 

‘  dA  j  dAj 

which  is  a  statement  of  Casligliano's  Theorem,  Part  I. 

For  purposes  of  later  developments,  it  should  be  emphasized  that  the 
strain  energy  in  Equation  III-3,  although  given  in  terms  of  displacements ,  could  have 
originally  been  expressed  in  terms  of  the  stresses  and  strains  in  the  clement, 
which  arc  in  turn  a  function  of  the  displacements. 

The  condition  employed  in  the  derivation  of  Equation  I1I-3  --  restraint  of 
all  displacement  components  except  the  one  of  interest  --  is  the  condition  associ¬ 
ated  with  the  definition  of  a  stiffness  coefficient.  Thus,  if  the  strain  energy  of  a 
discrete  element  can  be  expressed  in  terms  of  the  element  node  point  displace¬ 
ments,  application  of  Equation  III— 3  would  result  in  the  direct  determination  of  the 
element  stiffnesses. 

In  general,  it  is  not  possible  or  convenient  to  determine  the  elastic 
strain  energy  explicity  in  terms  of  the  node  point  displacements,  and  the  approach 
taken  is  as  follows. 

An  assumed  displacement  component  for  the  element  can  be  written  in  the  form 

A  =  ai  ft  (x,  y,  ■/.)  i  a2f2  (x,  y,  z)  . aNfN  (x,  y,  z)  (III-4) 

Here,  attention  is  restricted  to  developments  where  the  number  (N)  of  undetermined 
coefficients  alsa2  .  .  .  .  ajj,  is  equal  to  the  total  number  of  node  point  degrees  of 

freedom.  The  expressions  f^  (x,  y,  z) ,  f 2  (x,  y,  z)  may  be  polynomials  (e.g.,  xy2, 
y^,  etc.),  trigonometric  functions  (e.g.,  sinx,  cos  y,  etc.)  or  may  take  other  forms. 

Once  the  displacement  functions  are  chosen,  or  established  as  will  be  the 
case  when  stresses  or  strains  are  first  assumed  and  then  integrated  to  obtain  the 
corresponding  displacements,  they  can  be  evaluated  at  each  node  point,  resulting  in 
algebraic  relationships  between  the  node  point  displacements  and  the  constants 
ai,  a2,  .  .  .  a The  complete  set  of  such  relationships  for  all  node  points  of  an 
element  can  be  expressed  as 


(III-l) 

-1  can  be  written  as: 

(III—  2) 

(III-3) 
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W  MM 


(III- 5) 


who  re 


is  a  column  matrix  of  N  node  point  displacements. 

is  a  column  matrix  of  the  constants  aj_,  .  .  .  ajj. 

the  square  matrix  of  coefficients  in  the  relationships  between 

|  A  |  and  |  a  j  . 


Solving  Equation  III-5  through  inversion  of  f  B  j  ,  one  obtains 
[»}  [“!  "‘{A} 

Since  relationships  between  the  element  node  point  forces  ( 


(III- 6) 


the  displacements 


are  desired,  the  remaining  steps  in  this  development  per 


tain  to  the  establishment  of 
pose  consider  Equation  III-; 
written  as 


if  relationships  between  j  F  j  and  ja  j  .  For  this  pur- 
-3  which,  for  a  particular  node  point  force  Fj,  can  be 


dV_  d_F  oal  (  _dU_  ua2  ^  du  od  N 

i  d  A .  da,  dA.  da„  dA  .  '  da  dA. 

l  1  l  2  i  N  i 

(HI-  7) 

Application  of  III—  7  to  all  node  point  forces  results  in  N  relationships;  these  can  be 
summarized  in  matrix  form  as 

W  “  [^7]  {1^}  <ln-8' 


and,  as  shown  in  Reference  5  , 


[7^]  ' 


(III- 9) 


Also,  the  strain  energy  (TJ)  is  a  quadratic  form  in  the  stresses,  strains,  or  displace¬ 
ment  derivatives  and  is  therefore  a  quadratic  form  in  the  displacements  aj_,  ■  ■  ■ 


g  (a.  x  a.) 


(Ill- 10) 


Thus,  when  the  strain  energy  is  differentiated  with  respect  to  a  particular  constant, 
say  aj,  the  resulting  expression  is  a  linear  function  of  the  constants  a^,  .  .  an,  .  . 
a^,  and  can  be  written  in  matrix  form  as 
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-T^  *  l  °i  J  {  a  }  «  ■  l'  '  '  •  N> 

and  for  the  complete  set  of  derivatives  of  U,  one  obtains 

{tnr}  '  [c]  {“}  <m-12' 


Now,  by  combination  of  Equations  II-6,  -8,  -9,  and  -12,  there  results 


{F} 

■  u»nTM 

H  W 

(III— 13) 

w 

M  W 

(III— 14) 

M 

([b]-‘)t[c 

]  H'1 

(III- 15) 

where,  in  recapitulation 

£  B  "j  =  a  matrix  relating  the  node  point  displacements  to  the  undetermined 
constants  of  the  assumed  displacement  functions. 

C  j  =  a  square  matrix,  each  row  of  which  contains  the  coefficients  of  the 
constants  a^,  .  .  .  a^  in  equations  which  represent  derivatives  of  the 
strain  energy  with  respect  to  a^,  .  .  .  ajvj. 

Note  that  the  subscript  E,  used  in  the  previous  chapter  to  designate 
element  stiffness,  displacements,  etc.,  is  discarded  in  the  present  chapter. 


2.  Derivation  ol  Terms  Representing  Initial  Strain  Effects 

Initial  strains  can  be  defined,  for  present  purposes,  as  strains  that  exist 
in  a  structure  prior  to  the  imposition  of  applied  loads.  Only  initial  strain  due  to 
temperature  change  and  prior  inelastic  deformation  are  of  interest  to  this  report. 


The  total  value  of  a  strain  component  at  a  point  can  be  written  in  the  form 


e 


T  “ 


where  <r.p  is  the  total  strain,  is  the  strain  due  to  stress,  and 

initial  strain.  For  temperature  change 


(III— 16) 
e  1  is  the 


l 

e  = 


=  aT 


(III-17) 


and  for  accumulated  inelastic  strain 


i 

e  =  e 

P 


(III— 18) 
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The  development  of  a  procedure  for  deriving  terms  representing  initial 
strain  effects  requires  a  re-examination  of  the  expressions  for  strain  energy.  The 
strain  energy  to  be  employed  is  the  elastic  strain  energy,  but  the  strains  appearing 
in  the  related  formulation  must  be  total  strains  since  it  is  these  which  correspond 
to  the  displacements  of  the  element  stiffness  equations. 

Thus,  the  strain  energy,  U,  can  be  written  as 

i  rv 

U  =  -  /  fe  crd  V  (III—  1 9) 

and  by  substitution  of  the  expression  for  elastic  strain  (Equation  III—  1 G) 

»V 

U  |  j  (  C  T  -  c  *)  cr  d  V  (III-  20  ) 


The  significance  of  these  factors  will  be  shown  in  the  development  to  follow. 

If  the  formulation  of  element  properties  is  to  be  based  on  assumed  dis¬ 
placements,  the  relationship  between  the  undetermined  constants  and  the  displace¬ 
ments  is  simply  that  given  by  Equations  III — I  and  III— 5 .  The  corresponding  total 
strains  are  derived  by  differentiation  of  the  displacement  expression  and  the 
stresses  are  obtained  by  use  of  Hooke's  Law  (Equation  III-lGa).  From  equation  III— 20 , 
the  strain  energy  is 

E-  rv  i  2 

U  TT-  (  fT-  *)  dv  (III- 21  > 


or,  in  expanded  form 


(III— 2 2  ) 


and,  as  was  noted  earlier,  this  leads  to  a  quadratic  expression  in  the  a's  when 
the  assumed  displacement  function  is  substituted  and  the  integration  performed 
i.e., 


g(a^  x  a^)  +  h(a^  x  e  )  +  j  (  n}— 


(III—  23  1 


From  Equation  III—  8  and  III  -9 


{-,}=  (WV  {-%} 


(III- 24) 


£  b]  is  unchanged  by  the  presence  of  initial  strain  but  U  is  now  given  by  Equation 
(III- 10a),  rather  than  by  Equation  (III- 10).  Therefore 
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3.  Incremental  Stiffness 

When  midplane  stresses  or  forces  influence  the  behavior  of  plates  and 
bcamr  in  bending,  their  effect  on  the  analytical  formulation  of  flexural  stiffness 
relationships  is  in  the  form  of  an  "incremental”  stiffness,  i.e.,  as  an  addition  to  the 
usual  flexural  stiffness.  The  purpose  of  this  section  is  to  extend  the  proceeding 
formulations  to  include  techniques  for  incremental  stiffness  derivation.  This  develop¬ 
ment  applies  only  to  beams  and  plates  in  flexure  which  are  subjected  to  a  previously 
applied  and  equilibrated  known  axial  or  midplane  force  system. 

For  beams  and  plates  under  the  conditions  of  interest,  the  following 
relationships  exist  between  the  work  (W)  done  by  the  lateral  and  midplane  loads 
during  bending  deformation  and  the  strain  energy  of  bending  (Uj.)  (see  Reference  10), 

W  =  Wy  +  Wh  =  U  (III— 3  0) 

where 

=  work  done  by  the  lateral  loads  (F^.F^  .  .  .  F.  ...  F^)  during  flexure. 

=  work  done  by  the  midplane  loads  during  the  displacement  of  the 
structure  caused  by  the  lateral  loads. 

Furthermore,  the  change  in  W  with  respect  to  the  displacement  Aj  can  be  obtained 
directly  from  Equation  (III- 30)  as 
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aw 

v 


a(uf-wh 


(III- 31) 


3A;  aA. 

1  1 

Now,  if  the  restraint  condition  employed  in  the  formulation  of  Equation  III— 3  is  con¬ 
sidered,  it  follows  that 

W  VV 

v 


and  aw  dWv 

a a .  3a . 

I  1 

Consequently,  from  Equations  III- 3 ,  -31  and  -32 


F.  = 
i 


<Uf  -  wh) 


dA 


i 


(III— 3  2) 


(III— 33) 


The  transformation  of  Equation  III-33  into  a  matrix  formulation  of  the 
desired  element  stiffness  matrix  could  be  accomplished  rigorously  through  applica¬ 
tion  of  the  same  procedures  that  led  to  Equation  HI-15.  For  brevity,  however,  the 
already-developed  Equation  III— 15  will  be  used  as  a  basis.  First,  it  must  be  noted 
that  Wft  can  be  expressed  in  terms  of  the  lateral  displacements.  This  is  also  the 
form  taken  by  U^.  Hence,  (U^-W^)  can  be  defined  as  an  "effective"  strain  energy, 

U'  ,  and  Equation  III— 3  can  be  written  as 


d  U1 


(III— 34) 


Equation  III-34  is  of  similar  appearance  to  Equation  ni-3,  which  provided 
the  basis  for  the  development  of  Equation  III- 15.  To  use  Equation  111-34  in  the  same 
way,  one  must  note  that  of  the  two  matrices  making  Up  Equation  III-15,  only  the 
matrix  [c]  pertains  to  strain  energy.  Each  row  of  [c]  represents  an  equation  for 
the  derivative  of  the  strain  energy.  Since  now  the  strain  D  1  is  composed  of  two  parts 
one  can  write 

[°]  =  [Cf]  -[°h]  <in-35> 

in  which  ^  C^J  and  [c^j  result  from  the  required  operations  on  Uj  and  W,  respectively. 
Equation  III- 15  therefore  becomes 

M-N-W-fHTNM'1  -(t»nT[cJ[-r  <>«-> 
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It  is  pertinent  to  note  that  i'or  beam  segments 
F  /o  L  2 

wh  =  --r  L  m  -» 


(III— 3  7) 


and,  for  plates  of  constant  thickness 

wh [N*(ir)  1  Mif)2*2M^)(|R]dA<ra-38> 


4.  Stress-Displacement  Relationships 

It  is  possible,  in  theory,  to  determine  the  stresses  in  a  matrix  dis¬ 
placement  analysis  by  utilizing  the  solved-for  displacements  and  the  element  stiff¬ 
ness  matrices  to  evaluate  the  node  point  forces.  Then,  these  forces  are  transformed 
into  stresses.  It  is  believed  more  convenient,  however,  to  form  directly  a  set  of  re¬ 
lationships  between  the  element  stresses  and  node  point  displacements.  Such 
relationships  are  termed  "stress-displacement”  equations  in  this  report  and  are 
written  in  the  form 


W=  [«]W 


(III— 39) 


The  procedux-es  used  to  determine  these  equations  take  one  of  two  forms , 
dependent  on  whether  the  derivation  of  the  element  stiffness  properties  was  based  on 
assumed  stress  or  assumed  displacement  behavior.  Consider  first  the  case  of 
assumed  stress  patterns.  Here,  the  development  of  element  properties  begins  with 
expressions  of  the  form 

{^}=[D]  {a}  (HI- 40) 

Eliminating  the  column  of  constants  from  III-40  by  use  of  III— 24  results  in 

{  „  }.  [d]  [b]  -1  {a}  -  [d]  [b]-1  {  A1}  (111-41) 

or 

{  ^  }  =  [s]  {a}  -{ct1}  (HI— 42) 

where 

[s]  =  [d]  [b]  (HI— 43) 

{  <j  1 }  =  [  D]  [b]  _1  {  A1}  (HI— 44) 

The  initial  stresses  {  ct-1}  correspond  to  the  stresses  that  occur  when  the  node 
point  displacements  are  zero. 

When  the  derivation  of  element  properties  is  based  on  assumed  displace¬ 
ment  behavior,  the  starting  point  is 
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{a}  [B]  '»  {a} 


(III— 6) 


The  strains  must  be  derived  from  the  displacements  through  application  of  the 
strain-displacement  derivatives.  When  this  differentiation  is  applied  to  the  basic 
assumptions,  one  obtains 


{•}  M  W 


and,  from  Hooke's  Law 


W  M  {«} 


Hence,  by  combination  of  III— 6,  -  45  and  -  46 


{'}-  M  MM'lM  -  {*'} 


thus,  in  the  present  case 


[s]  [K]  [H]  [b]  -1 


(III-45) 


(III— 46) 


(III-47) 


(III-48) 


B.  BEAM-AXIAL  FORCE  ELEMENT 


1.  Basic  Considerations 

The  beam-axial  force  element,  shown  in  Figures  III-2  and  III— 3 ,  appears 
in  each  of  the  three  computer  programs  described  in  Chapters  IV- VI.  It  is  the  only 
major  element  in  the  One-Dimensional  Structure  Program  (Chapter  IV),  where 
detailed  attention  is  given  to  the  variation  of  temperature,  etc.,  on  the  cross-section 
of  the  element;  only  the  gross  effects  on  element  cross-sections  are  treated  in  the 
other  two  programs.  In  developing  the  relationships  for  this  element,  the  expres¬ 
sions  employed  in  the  One- Dimensional  Structure  Program  will  be  established. 

By  introducing  the  simplifications  pertinent  to  the  other  two  programs,  the 
derived  relationships  can  be  reduced  to  the  forms  of  interest. 

This-dey-elopment  is  based  on  the  following  assumptions: 

(1)  Cross-sections  originally  plane  remain  plane. 

(2)  The  geometry  and  temperatures  on  a  cross-section  are  symmetric 
about  one  axis;  the  bending  moments  act  only  in  the  plane  of  symmetry. 

(3)  The  cross-sectional  geometry  and  temperatures  do  not  vary  along  the 
axis  of  the  element. 

(4)  Plastic  strains  are  constant  over  the  length  of  the  element  and  are 
dependent  only  on  the  stress  conditions  on  a  cross-section  midway 
between  the  ends  of  the  element. 

(5)  Lateral  deflections  are  relatively  small. 
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y 
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Figure  iri-2.  One- Dimensional  Element 


Figure  III— 3.  Typical  Cross-Section  of  One-Dimensional  Element 
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Consider  a  typical  symmetrical  cross-section  of  the  eiemont  as  shown  in 
Figure  III— 3.  In  accordance  with  assumption  (1),  the  total  axial  displacement  (u,p)  of 
an  arbitrary  point  on  the  cross-section  can  be  written  as 


U T  UT  ‘  9T  t 


(III-49) 


where  £  is  the  distance  to  the  point  in  the  direction  of  the  axis  of  symmetry  measur  ’d 
from  an  arbitrarily  chosen  reference  lint  and  is  both  the  angular  displacement 
of  the  cross  section  and  the  slope  of  each  fiber  of  the  beam  at  the  cross-section  under 
study.  "Barred"  values  for  u  and  e  represent  quantities  associated  with  pure  trans¬ 
lation  of  the  cross-section.  By  definition: 


9. 


dw.j, 

dx 


Since  the  total  axial  strain,  e  x 


'T  dx 


9 

d“\V-r 


dx  2 


(III— 50) 


T 


^T 

dx 


XT 


1 
P  x 


-  ' '-fill-  51) 


The  term  £  /  p  ,p  represents  the  component  of  the  total  strain  re  salting  from  the 
relative  rotation  of  the  cross-section,  where  P  x  is  the  total  cu?vature. 

The  total  strain,  e  x  is  composed  of  two  components:  (1)  an  elastic- 

strain  — and  (2)  an  initial  strain  €  Replacement  of  the  total  strain  in 
F, 

Equation  III-49  bv  the  sum  of  its  components  results  in 


(111-51) 


It  is  now  necessai'y  to  relate  the  stresses  to  forces  and  moments  cr 
’Jstress- resultants".  The  axial  stress  resultant,  F  ,  and  the  bending  stress  resultant, 
My,  are  related  to  the  stresses  on  the  cross-section  through  equilibrium  as  follows 


F 

x 


M 

,y 


(III- 53) 
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(The  "bar" 

over  My  signifies  an 

internal  moment. 

See  Chap.  11), 

and,  by  use 

of  Equation  III-52. 

F  = 

X 

e  f  EcLA  - 

*T  ' 

/•A 

/  h£  dA  - 

PT 

rA 

J  E  e  1  c!A 

J  X 

M 

y 

<F  f  E  ^  tIA  - 

XT 

fA  E  C  2  ,1A  / 

PT  J 

■A 

£  E  c  1  dA 

X 

(III- 53a) 

These  equations,  can  be  simplified  by  referring  the  integrations  to  the  elastic  neutral 
axis  location ,  £  ,  defined  as 


£ 


fA  AK  £  dA 
fA  E  dA 


(III—  54) 


Hence,  if  the  variable  £  in  Equations  (III-53a)  is  replaced  by  £  £  -  £  (i.e.,  the 

origin  of  coordinates  is  now  placed  at  the  neutral  axis),  J  A  E  £  '  dA  =  0  and  the 
equations  reduce  to 


F  e 

x 


EA  -  F 


xT 


M 


-  M  1 

p  T  y 


(III—  53b) 


or 


F  +  F 
x  x 

EA 


F 

x 

EA 


EA 


ex 


EA 


(M  *  M  \ 

_ XJ 

El 


M 

_X 

El 


M  1 


El 


1 

Px 


M 

1j 


(ITI— 53c) 


whe  re 


/-A 

EA  = 

J  EdA  (elastic  axial  stiffness) 

El  =  , 

r* 

(£’  >2 

dA  (elastic  flexural  stiffness) 

rA 

F  1  = 

X 

J  E 

dA  (initital  axial  force) 

M*  =  . 

rA° 

dA  (initital  moment) 

(111-55) 
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e  Average  elastic  strain  due  to  the  axial  stress  resultant, 
x 

e 

P  e  curvature  due  to  the  bending  stress  resultant. 

Thus,  the  integrals  containing  the  initial  strain  terms  are  in  the  form  of  initial 
equivalent  stress  resultants.  Note  that  in  replacing  £  by  £'  ,  the  displacement 
u  is  now  defined  as  being  the  displacement  at  the  neutral  axis.  Also, 

i  dil 

€  C  _ T 

XT  XT  dx 
and,  (from  111—52) 


Use  of  Equations  (III-53c)  in  Equation  (111-52)  results  in 


cr  E 
x 


(F  ■>  F  S  (M  '  M  ‘)  e 

— 5 - S-  +  - 1 — —X — *  _  e 

EA  El  x 


(III-52a) 


Equations  III-52a  and  111-53  are  the  basic  equations  for  this  development. 

It  is  important  to  note  that  the  location  of  the  neutral  axis  as  defined  by  Equation 
III-54  is  independent  of  the  magnitude  or  distribution  of  the  initial  strains;  it  depends 
only  on  the  variation  of  Young's  Modulus  (E)  which  in  turn  is  only  a  function  of 
temperature.  Thus,  the  axial  stiffness  EA  and  the  flexural  stiffness  El  only  depend 
on  the  cross-sectional  geometry  and  the  temperature  and  if  E  is  a  constant  the 
neutral  axis  will  pass  through  the  centroid  of  the  cross-section.  It  should  also  be 
noted  that  the  axial  force  Fx  is  an  axial-stress  resultant  with  its  point  of  applica¬ 
tion  at  the  neutral  axis  on  the  axis  of  symmetry.  Correspondingly,  the  bending- 
stress  resultant  (My)  acts  about  the  neutral  axis.  Thus,  externally  applied  axial 
forces  are  presumed  to  act  through  the  neutral  axis  of  the  cross-section. 

2.  Axial  Behavior 

Since  the  conditions  associated  with  axial  load  (see  Figure  III—  2)  pro¬ 
duce  a  state  of  constant  strain  in  the  axial  direction,  the  axial  displacements  are  a 
linear  function  of  the  axial  coordinate  and  points  on  the  neutral  axis  can  be 
expressed  as 


+  a^x 
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w 


From  Equation  (III-12a),  it  follows  that 


(111-61) 


and  from  Equation  (111-26)  through  (111-29) 


3.  Force-Displacement  Equations--Flexural  Behavior 


In  the  ease  of  flexure  (see  Figures  I II— 2  and  III-4)  a  distinction  must  be 
made  between  the  initial  displacements  associated  with  thermal  and  plastic  strain 
and  the  initial  displacements  produced  during  fabrication.  This  is  necessary  be¬ 
cause  the  latter  are  exempt  from  geometric  boundary  conditions,  i.e.,  if  it  is 
specified  that  a  node  point  displacement  lx.'  zero,  it  is  meant  that  only  the  displace¬ 
ment  due  to  applied  loads,  temperatures  and  plastic  strain  is  zero — the  node  point 
displacement  due  to  fabrication  remains.  To  retain  this  distinction,  the  total 
transverse  displacement  (w,,,1)  is  considered  to  be  composed  of  three  parts: 

w  .j,1  w  '  w.  •  w  (II1-G4) 

who  re 

w  the  transverse  displacement  due  to  elastic  strain 
e 

w  .  the  transverse  displacement  due  to  thermal  and  plastic  strain 

w  f  the  fabricational  transverse  displacement 


Figui'e  III— 4 .  Beam  Segment  with  Fabricational  Displacements 
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Generally  w'-p  (the  total  displacement  from  a  reference  axis)  is  employed 
in  problems  involving  initial  fabricational  displacement.  Here,  however,  it  has 
been  decided  to  employ  the  displacements  corresponding  to  the  difference  Ixdween 
the  total  and  fabricational  displacements  in  the  definition  of  the  force-displace¬ 
ments  relationships  in  the  program.  These  displacements. w-p  in  this  case,  are 
therefore  defined  as 


wT  (w'T  -  w  f)  (we  t  w.) 


(Ill—  65 ) 


To  derive  the  equations  of  interest,  the  displacements  are  assumed  to 
have  the  form 


WT  "I  ’  V  '  a3  X  "  a4X 


(III—  G  6) 


so  that 


dw,r  2 

B  n  •<  2a  x  +  2a  x 

T  dx  2  3  <1 


(111-67) 


T 


d2\v 

“2 


1  2a3  *  (ia4x 


(III- 68) 


By  evaluation  of  III-66  and  -67  at  the  node  points  there  is  obtained 
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>  (III- 69) 


V 


a 

4  J 


Consider  first  the  component  associated  with  simple  flexure,  given  in 
this  case  by 

•L  _ 

(in- 70) 


f  -  1 

J  M  (  —  )  dx 

o  y  p 


_ _  _  —  i  —  &  — —  n 

By  use  of  Equations  III-53b  and  III-53C,  and  with  M  M  +  M  this 

becomes  ■'  ' 
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From  Equations  (III-36)  and  (III-29),  and  with  use  of  (111-69),  it  follows  that 


N  -  <  H  ">T[ 

>.]  H  -1 

(III-36a) 

r  _  a )  ,  r  _  1  -l,  T 

r  dh  -|a 

{f,  }=  -  <  L  B  j  ) 

(III—  29a.) 

r  i  -i  t  ( 

■dh\  P 

n 

i 

X 

N 

.  da  } 

(Ill— 29b) 

The  resulting  matrices  (  [kf]  ,  { Fz  }  ,  {f  Pj)  are  shown  in  Figure  II1-5,  where 

the  subscript  T  lias  been  dropped  from  w  and  Q  . 


Figure  III—  5 .  One  Dimensional  Element  —  Out-of-Planc 
Force-Displacement  Relationships 


The  terms  of  the  force-displacement  relationships  that  arise  from  the  work 
done  by  the  axial  force  during  flexure  will  now  be  developed.  Here,  the  term  repre¬ 
senting  work  done  by  inplane  forces  is  (see  section  III. A. 3) 
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(III— 37a) 


where  Fx  is  the  axial  stress  resultant,  which  is  a  known  value  prior  to  the  flexural 
analysis.  Note  that  the  slope  due  to  the  total  displacement  from  the  reference  axis 
is  employed  in  the  above  equation  since  Equation  lll-37a  emanates  from  strictly 
geometric  considerations.  This  total  displacement  is  assumed  to  be  of  the  form: 

w'T  =  a,1  +  a'2  x  +  a'3  x2  +  a'4  x3  (III-72) 

Together  with  Equation  (111-6(1),  this  assumption  implies  that  the  initial  fabricational 
displacements  take  the  same  form  as  the  elastic  displacements. 

With  use  of  the  derivative  of  Equation  (III-72)  in  (lll-37a),  there  is  obtained: 

(aV2L+T(aV2L3+l(a'4)2L5  +  2a,2a’3L2 

(III-37b) 

+  2a’  a’  L3  4  3a’  a'  L,4  1 
2  4  3  4  J 

and,  by  establishing  the  derivatives  of  with  respect  to  each  of  the  constants 
a',  a1  there  is  obtained. 


F  f 

- _ x_ 

2 


Since  the  stiffness  relationships  are  given  by  (  TbI  1)T  [  C.  ]  (^bJ  1  this  parti¬ 
cular  contribution  to  flexural  stiffness  (denoted  by  [n  1  )  is  given  by 
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[„]=<[B]-V[ch]  [b]"1 


(III-74) 


The  postmultiplier  of  the  above  incremental  stiffness  matrix  is  the  column 
of  displacements  W  and  Q- p'.  As  stated  earlier,  however,  it  is  desired  that  all 
relationships  be  expressed  in  terms  of  the  difference  between  w'-p  and  the  fabri- 
cational  displacements  wfl  i.e.,  in  terms  of  wT^.  This  is  easily  accomplished  by 
making  the  substitutions  w'T  =  wT  +  w{,  Q  =  9  T  +  9  resulting  in 


dr 


9^ 


9  r 


T„ 


[»K 


f  \ 
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9  r 


>  '  [■] 
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T 


(III-74a) 


also,  since  the  fabrieational  displacements  Q  w  ^  are  known  quantities  in  an 
analysis,  the  product  of  the  incremental  stiffness  and  the  column  of  fabrieational 
displacements  (the  second  term  on  the  right  side  of  Equation  III-74a)  can  be  ex¬ 
pressed  as  an  initital  force  column  |  F^  j  ,  i.e. 


K'}  -  [»]J 


e, 


(III- 7 5) 


w  f 
V  2  ) 


The  explicit  forms  for  [  n]and  |  Fz  }  appear  in  Figure  (III- 5) ,  where  they  are 
included  in  the  complete  expression  for  the  out-of-plane  force-displacement 
equations  for  the  one-dimensional  element. 

4.  Stress- Displacement  Equations 


In  general,  the  stresses  will  vary  along  the  length  of  the  element  so  that 
a  presentation  of  the  complete  stress  distribution  would  involve  a  voluminous  amount 
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of  information.  Thus,  only  the  stresses  at  the  center  of  the  element,  which  should 
suffice  to  define  the  significant  stresses  within  the  element  and  structure  as  a  whole, 
are  determined. 


Equation  (III-52a)  gives  the  stress  at  a  point  on  the  cross-section  of  the 
one-dimensional  element  in  terms  of  the  stress  resultants.  The  stress-resultants 
have  been  formulated  in  terms  of  the  displacements  in  the  preceding  two  sections. 
By  combining  these  formulations  with  Equation  (III-52a)  and  effecting  an  evaluation 
at  the  central  cross-section  of  the  element  (i.e.,  at  x  L/2)  the  following  stress- 
displacement  equation  is  obtained 


<T 


X(x  L/2) 


Sy2>-<°  AT"  «/>] 
(111-76) 


C.  PLATE  ELEMENTS 
1.  Triangular  Plate 

A  detailed  development  of  all  plate  element  relationships  w'ould  be  beyond 
the  scope  of  this  report.  Their  development  is  documented  elsewhere  (Refei'ence  8). 
The  present  report  simply  outlines  the  bases  for  the  derivation  of  the  pertinent  stiff¬ 
ness  matrices,  thermal  forces,  etc. 

The  stiffness  matrix  for  the  triangular  plate  element  was  originally 
derived  by  Turner,  ct  al  (Refei'ence  3).  By  use  of  Castigliano's  Theorem  (Equation 
III- 15)  and  the  assumption  oT  constant  strains  8 x,  8 y,  Y\  the  stiffness  relation¬ 
ships  were  re-derived  for  orthotropic  behavior.  The  element  is  shown  in  Figure 
III— 6 .  Note  that  the  element  is  arbitrarily  oriented  in  the  x-y  plane;  use  is  not  made 
of  one  of  the  element  edges  as  a  local  coordinate  axis.  The  x-y  axes  are  intended 
to  be  the  axes  of  the  complete  structure.  This  is  because  it  is  expected  that  the  axes 
of  the  complete  structure  will  correspond  to  the  principal  axes  of  orthotropy.  The 
inplane  thermal  and  plastic  forces  were  derived  from  Equation  111-29  under  the 
assumption  that  the  thermal  and  plastic  strains  were  constant  thr-oughout  the  element. 

After  extensive  investigation,  it  was  decided  to  base  the  development  of  the 
relationships  for  out-of-planc  behavior  on  the  following  assumed  displacement  func¬ 
tion. 


2  2  3  3  2 

a  +a  xh  a„y  +a  y  >  ax  +  a„x  +  a  y  +a  xy  +  a  xy 
12  :r  4  5  6  T  8  9 


(111-77) 


It  is  to  be  noted  that  this  function  is  geometrically  unsymmetric,  i.e.,  the  term 
a jX^y  which  is  the  counterpart  of  a^xy2  is  absent.  Attempts  were  made  to  remedy 

this  situation  but  none  succeeded.  By  operating  on  Equation  (III-77)  in  the  manner 
indicated  by  Equation  III-36,  the  desired  flexural  and  incremental  stiffnesses  were 
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derived.  The  thermal  and  plastic  moment  terms  were  obtained  by  means  of  a 
prorating  of  the  "edge  moments"  to  the  corner  points.  Both  the  thermal  and  plastic 
moments  were  assumed  constant  throughout  the  plate. 

2.  Quadrilateral  Plate 

The  quadrilateral  plate  clement  is  shown  in  Figure  II1-7 .  The  relation¬ 
ships  for  plane  stress  behavior  arc  those  adopted  by  Turner,  et  al,  in  Reference 
3,  i.e. 


a  - 

X 

ai  +  V 

V 

a3  +  ‘V 

Txy  = 

a5 

(III- 78) 

There  are,  in  fact,  numerous  alternatives  to  these  assumptions.  In  the  extensive 
evaluations  of  Reference  5,  however,  it  has  been  found  that  there  are  unimportant 
differences  in  the  results  for  these  alternatives  when  other  than  a  coarse  gridwork 
of  points  is  involved.  By  use  of  the  above  assumptions  and  Equation  III- 1 5 ,  stiffness 
equations  were  formulated  for  the  case  of  orthotropic  material  properties.  With 
respect  to  thermal  and  plastic  forces,  it  was  decided  to  formulate  the  required 
expressions  for  condition  of  constant  thermal  and  plastic  sti'ain. 

A  suitable  basis  for  the  derivation  of  the  flexure  stiffnesses  was  not  to  be 
found  elsewhere.  The  approach  adopted  was  first  to  assume  the  following  displace¬ 
ment  function: 


w 


3  2  3  2  3 

=  a^x  +  a^x  +a^x  +  a^y  +  a y  +  a^y  +  a^x  y 

2  3  2 


+  agx  y  +  a9xy,  a1()xy 


ailXy 


12 


(III— 7  9) 


Then,  by  means  of  Equation  (III- 36)  the  appropriate  flexural  and  incremental  stiff¬ 
nesses  were  derived.  Note  that  (111-79)  is  "geometrically  symmetric",  e.g.,  cor¬ 
responding  to  the  agx^y  term  there  is  an  a-Q.xy^  term.  Also,  (III- 79)  satisfies  the 
appropriate  differential  equation  of  equilibrium.  To  obtain  thermal  moments,  the 
temperature  gradient  across  the  thickness  was  assumed  constant  at  all  points  on  the 
surface  and  corner  thermal  moments  were  defined  by  prorating  the  distributed  edge 
moments  to  the  comers. 
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Figure  III— G.  Triangular  Plate  Element 


v. 


Figure  III—  7.  Quadrilateral  Plate  Element 
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CHAPTER  IV 

PROGRAM  FOR  THE  ANALYSIS  OF  ONE-DIMENSIONAL 
STRUCTURES 

A.  SCOPE 

The  computer  program  described  in  this  chapter  provides  a  means  for  the 
determination  of  the  stresses,  displacements,  and  instability  of  structures  posses¬ 
sing,  cross-sectional  dimensions  which  are  small  with  respect  to  the  length  dimen¬ 
sion  and  where  the  deformational  behavior  on  all  cross-sections  is  governed  by 
elementary  beam  theory,  i.e.,  the  assumption  that  plane  sections  remain  plane  under 
deformation.  In  addition  to  applied  loads  and  thermal  gradients,  the  effects  of  initial 
displacements  and  plasticity  can  be  taken  into  account. 

Figure  IV-1  illustrates  typical  conditions  which  can  be  treated  with  use  of  the 
program.  The  "real"  structure  appears  in  Figure  IV-la,  while  the  idealization 
appears  beneath  in  Figure  IV-lb.  Each  discrete  element  is  a  one-dimensional 
segment  of  constant  cross-section  (Figure  lV-lc).  Comparison  of  Figures  IV-la  and 
IV-lb  discloses  the  limitations  and  capabilities  of  this  type  of  idealization.  Length¬ 
wise  variations  of  cross-sectional  geometry,  as  well  as  temperature  conditions,  are 
represented  in  a  stepwise  manner.  Distributed  loads  must  be  defined  in  terms  of 
statically  equivalent  concentrated  forces  at  the  element  juncture  points.  Support 
conditions,  including  flexible  support,  can  be  imposed  at  each  node  point.  The  repre¬ 
sentation  of  a  flexible  support  consists  of  "flexible  support"  (or  "restraint")  ele¬ 
ments.  To  distinguish  between  such  elements  and  usual  elements  (pictured  in  Figure 
IV-lc)  whenever  confusion  is  likely  to  result,  the  latter  will  be  referred  to  as 
"conventional"  elements. 

The  force-displacement  behavior  of  a  conventional  element  is  segregated  into 
inplane  and  out-of-plane  behavior,  respectively.  Formulation  of  the  detailed  rela¬ 
tionships  for  this  element  was  accomplished  in  the  previous  chapter. 

The  maximum  of  fourteen  conventional  elements  and  four  flexible  support  ele¬ 
ments  can  be  employed  in  the  idealization  of  a  structure.  Initial  displacements  of 
the  structure  due  to  fabrication  are  taken  into  account  by  specifying  the  values  of 
these  displacements  at  the  element  juncture  points. 

In  order  to  provide  for  an  arbitrary  (but  symmetric)  geometric  and  tempera¬ 
ture  condition  on  the  cross-section  of  each  element,  provision  is  made  for  the  cross- 
section  of  each  element  to  be  divided  into  as  many  as  30  differential  areas  d  A. 

(See  Section  A-A  of  Figure  IV-lc).  A  value  of  temperature  can  be  assigned  to  each 
differential  area,  thereby  defining  the  cross-sectional  variation  of  temperature  and 
also  of  material  properties,  since  the  latter  can  be  defined  as  being  temperature 
dependent. 
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a.  Actual  Structure 


b.  Analytical  Idealization 


Z 


Axis  of 


c.  Individual  Element 


Figure  IV-1.  Conditions  for  One-Dimensional 
Analysis  Program 
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Flexible  support  elements  are  described  by  their  axial  stiffness  (  r  ),  lateral 

stiffness  (  £w).  and  flexural  (or  torsional)  stiffness  (  )•  These  quantities  will 

also  be  referred  to  as  "restraint  coefficients".  Each  flexible  support  element  is 
associated  with  two  node  points  and,  in  general,  displacements  of  the  node  point  of 
the  restraint  element  which  is  not  attached  to  the  structure  proper  is  considered  as 
fixed.  Support  conditions  at  the  other  (attached)  node  point  must  be  specified  in  con¬ 
formity  with  the  conditions  of  the  problem  at  hand. 

In  general,  the  input  must  consist  of  the  data  for  every  element,  including  all 
differential  areas  and  their  corresponding  temperatures.  However  there  are  many 
practical  circumstances  in  which  simplified  conditions  prevail  and  for  certain  of 
these  a  provision  has  been  made  for  reduced  input.  These  options  are  detailed  in 
the  report  describing  input  data  preparation  for  the  program. 

The  capabilities  for  inelastic  analysis  permit  solutions  for  either  time- 
independent  or  time-dependent  (creep)  behavior,  or  both  in  combination.  The  rela¬ 
tionships  for  such  analyses  are  discussed  in  the  next  section,  where  the  theoretical 
basis  of  the  program  is  outlined.  The  computational  process  for  inelastic  analysis, 
which  differs  considerably  from  the  process  for  elastic  and  instability  analysis,  is 
treated  in  a  later  separate  section,  however, 

B.  THEORETICAL  BASIS 


1.  Elastic  Analysis 

The  following  is  a  description  of  the  analysis  procedure  for  conditions  of 
linear  elastic  behavior,  including  thermal  stress,  initial  displacements,  and  elastic 
instability.  The  terms  for  inelastic  analysis  are  represented  since  the  procedure 
for  inelastic  analysis  is  essentially  a  succession  of  elastic  analyses.  (The  proce¬ 
dures  for  inelastic  analysis  are  examined  in  the  next  section). 

As  shown  in  the  preceding  chapter  (Section  IB.B)  the  force  displacement 
relationships  for  the  beam-axial  force  element  are 


For  inplane  behavior 


k,}  ■  kl  W  *  K)  *  k, 


n 


(iv- 1) 


For  out-of-plane  behavior 


{m  , 

F  ] 

r  k  i 
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z  J 
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l  ZJ 

*  K.  f,  }  ♦  { 
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Where 


{  Fx}  •  {My’FJ 

are  column  vectors  which  list  the  forces 
acting  upon  the  ends  of  the  element, 

[ 

are  the  respective  elastic  stiffness 
matrices 

{  ».  },  {  ey,w} 

are  column  vectors  which  list  the  dis¬ 
placements  at  the  ends  of  the  element. 

{  P‘°}  ■  { K  ■ Fz°  } 

are  forces  produced  by  complete  restraint 
of  the  element  against  thermal  deforma¬ 
tions. 

{  f/Hs/^p} 

are  forces  produced  by  complete  restraint 
of  the  element  against  inelastic  deforma¬ 
tions. 

{  S/’  F/} 

are  equivelent  forces  dependent  on  the 
presence  of  fabricational  displacements 
at  the  ends  of  the  element. 

[”*] 

"Incremental  Stiffness"  -  the  influence  of 
the  axial  loads  on  flexural  behavior. 

The  detailed  form  of  these  matrices  was  presented  in  Chapter  III.  There,  it  was 
shown  that  each  matrix  on  the  right  side  of  Equations  (IV-1)  and  (IV-2)  contains  a 
scalar  multiplier  which  is  an  integral  to  be  evaluated  over  the  area  of  the  cross 
section.  k  j  .  for  example,  contains  the  scalar  /  E*dA  ,  while.  £  k  J  has  the 

multinlier  2/£2E»dA 

Z3 

The  evaluation  of  the  above  integrals  is  accomplished  through  a  subdivi¬ 
sion  of  the  cross  section  of  each  element  into  differential  areas.  As  noted  earlier, 
the  element  cross  section  can  be  subdivided  into  as  many  as  thirty  differential 
areas;  the  temperature  variation  on  the  cross  section  is  represented  by  assigning 
an  average  temperature  to  these  areas.  Since  the  material  properties  are  repre¬ 
sentable  as  a  function  of  temperature,  the  computer  selects  for  each  differential 
area  those  properties  which  are  consistent  with  the  temperature  of  that  differential 
area. 


Based  on  geometric,  load,  temperature,  and  material  property  data  of 
the  problem  and  the  selected  analysis  options  the  computer  evaluates  the  pertinent 
portions  of  Equations  (IV-1)  and  (IV-2)  for  each  element.  The  individual  element 
relationships  are  then  assembled,  in  satisfaction  of  node  point  equilibrium  and  com¬ 
patibility  requirements  (See  Chapter  II),  and  geometric  boundary  conditions  are 
applied  to  yield  a  set  of  equations  for’  the  assembled  "analytical  model". 
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Once  the  stiffness  equations  have  been  formulated  the  analysis  process 
can  take  one  of  many  forms,  dependent  upon  the  exercise  of  available  options,  but  in 
the  following  discussion  it  will  be  assumed  that  none  of  the  possible  analysis 
operations,  except  those  pertaining  to  inelastic  behavior,  are  to  be  eliminated, 
(Analysis  for  inelastic  behavior  is  discussed  in  the  next  section).  In  such  a  case, 
the  relationships  for  the  complete  structure  can  be  written  in  matrix  form  as: 


For  axial  behavior 


(p  1  =  |~K  1  ful  +  (p  Q1  +(pPl 

lxJ  LXJIJ  lxJ  lxJ 


(IV— 3) 


For  flexural  behavior 


{*VP«}'K]  {VWHN]  h’"}  *{My  ,P*°}  <IV-4> 

*  w-p4 

These  equations  are  presented  in  general  form  in  Chapter  in  as  Equations  (11-12)  and 

(II  l3)-  In  fjrst  step  0f  solution  process  the  axial  displacements,  ju  j  , 
are  determined  by  inversion  of  the  r  j  matrix  in  Equation  (TV-3). 


Thus: 


{“}  ■  [Kx]_1{v  Px°-PxP} 


av-5) 


These  axial  displacements  ( j  u  j  )  are  substituted  back  into  the  element  relation¬ 
ships  (Equation  IV-1)  to  yield  the  element  node  point  forces  ( {  Fx  }  )>  from  which 
the  axial  stress  resultants  are  obtained. 

The  influence  of  the  axial  forces  is  accounted  for  in  the  flexural  analysis. 
The  matrix  Jn  j  is  a  function  of  the  axial  loads  and  can  be  evaluated  after  the 
inplane  analysis  is  complete.  The  solution  for  the  flexural  displacements  follows 


(VWi  “  [KJ  *Hl‘1{{My'PJ-{My°  'PxJ  «V-6> 

-{My’pPHMy'PI}} 

These  flexural  displacements  and  the  previously  determined  axial  stress  resultants 
are  employed  in  the  direct  determination  of  the  stresses  on  each  of  the  differential 
areas  of  each  structural  element. 

To  accomplish  a  determination  of  a  critical  load  it  is  first  necessary  to 
specify  a  Reference  Critical  Load  Value-Prejv  All  element  axial  loads  are  nor¬ 
malized  on  this  value  through  the  relation 
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X  P  . 

n  ref 


av-7) 


where  Fn  is  the  axial  stress  resultant  on  the  nth  structural  element. 

Equation  (IV-4)  then  appears  as: 

{ M  ,P  }  =  f  K  1  f  fl  ,w  U  P  Jn  1  (  9  ,w  )  +{m  °  ,  P 

I  v  zl  L  ■/.  1  l  v  J  rei  l  J  l  y  ■  J  ^  y  z  i 


+|mp 

y 


,  pp  1  h  {  Mf  ,  Pf  j- 
z  .  I  y  /.  J 


(1V-8) 


N  I  now  contains  the  X  values,  rather  than  the  F  values. 

[I  n  n 

The  condition  of  clastic  instability  derived  from  Equation  IV-8  is 


p~  tvw}  '  K1  ‘VI  tVwl 

ref  J 

Application  of  matrix  iteration  to  tins  condition  yields  the  minimum  P  ^  and 
corresponding  mode  shape  for  instability. 

With  regard  to  the  selection  of  the  reference  axial  load,  the  program 
allows  for  two  possibilities.  In  one,  the  reference  axial  load  is  automatically  set 
equal  to  the  axial  force  computed  in  a  designated  element.  Alternatively,  the  value 
of  the  reference  load  can  be  designated. 

Once  the  critical  reference  axial  load,  PreI-  lias  been  computed,  the 
distribution  of  critical  axial  forces,  as  obtained  from  Equation  (IV-7),  is  given  by: 


F 

n 


cr 


x 


n 


ref 


cr 


(IV-10) 


Neither  the  foregoing  discussion  nor  Chapter  111  has  given  explicit  con¬ 
sideration  to  the  treatment  of  flexible  support  elements.  Flexible  support  conditions 
at  a  node  point  are  simulated  by  means  of  this  special  type  of  element  which,  as  in 
the  case  of  the  conventional  element,  is  associated  with  two  node  points  and  three 
stiffnesses,  one  each  in  the  direction  of  the  three  displacement  degrees  of  freedom. 
These  three  stiffnesses  (or  "restraint  coefficients"),  £ £w,  and  4,  ,  are 
defined  as  follows: 

For  axial  behavior: 

AF  =  (C  )  Au  av-ll) 

x  u 
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For  lateral  behavior 

AFz=  (4W)  Aw  (IV- 12) 

and,  for  rotational  bchavior 

AM  =  (4.)Af?  (IV-13) 

y 

Note  that  the  three  restraint  coefficients  are  independent  of  each  other.  It  is  left  to 
the  analyst  to  decide  upon  the  values  of  these  coefficients.  By  analogy  with  the  stiff¬ 
ness  matrices  for  a  conventional  element,  the  stiffness  matrices  for  a  restraint 
element  are: 


For  axial  behavior: 

K1 

and,  for  flexure. 

[  -,] 


C,  -c. 


o  o 
o  o 


0 

0 

-  c, 


(IV-14) 


(IV-15) 


Effects  of  initial  displacements  due  to  fabricational  inaccuracies 


accounted  for  by  the  inclusion  of  a  column  matrix, 


{S'.w'} 


,  composed  of 


initial  rotational  (  Q  )  and  lateral  (w  )  deflections  determined  at  each  node  point. 
These  initial  deflections  are  measured  relative  to  a  straight  line  which  passes 
through  the  x-axis  of  the  one-dimensional  structure.  Total  deflections  are  deter¬ 
mined  as  the  addition  of  the  fabricational  deflections  and  the  deflections  computed 
from  Equation  IV-G. 


2.  Inelastic  Analysis 

It  is  intended  that  this  program  be  capable  of  solving  inelastic  problems 
wherein  the  loads  and  temperatures  vary  with  time,  as  sketched  below.  For  analyti¬ 
cal  purposes  and  to  avoid  the  necessity  of  listing  large  quantities  of  time  increments 
that  may  be  required  in  the  computations,  points  in  a  given  load-temperature  history 
are  first  specified.  These  points  are  illustrated  by  the  heavily  dashed  lines  in  the 
sketch  below;  the  distance  between  points  are  termed  time  intervals  in  this 
report.  Each  time  interval  may  then  be  further  subdivided  into  equal  time  incre¬ 
ments  which  may  differ  for  each  interval. 


In  an  earlier  development  of  the  matrix  displacement  method  to  permit 
inelastic  analyses  (Reference  9)  the  approach  suggested,  was  based  on  the 
selection  of  time  increments  small  enough  so  that  a  single  computation  cycle  for 
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(or  pertaining  to)  the  time  increment  would  suffice  to  determine  accurately  inelas¬ 
tic  strains,  both  time  independent  and  creep  strains.  However,  to  afford  greater 
flexibility  and  possibly  more  accuracy  under  certain  conditions,  it  is  proposed  here 
that  time-independent  plastic  strains  be  computed  at  the  ends  of  each  time  increment 
(i.e.,  at  specific  times  in  the  history)  by  either  an  accumulative  step-by-step  or  itera¬ 
tive  process.  Creep  strains,  on  the  other  hand,  are  computed  within  a  time  incre¬ 
ment  based  on  the  assumption  of  constant  stress  and  average  temperature  conditions. 
For  creep  determinations  the  stress  is  taken  equal  to  the  stress  which  prevails  at 
the  start  of  a  time  increment. 

There  are  no  limits  on  the  number  of  time  intervals  and  increments  which 
can  be  used  in  an  analysis;  the  number  of  accumulative  steps  or  iterations  performed 
at  a  given  time  in  the  load-temperature  history  in  an  inelastic  analysis  are  limited, 
however,  by  a  specified  convergence  criterion  with  regard  to  the  stresses. 


Load 


Time  (t) 

The  time-independent  plastic  strain  component  of  the  total  strain  €  ^  obtained 

in  a  conventional  uniaxial  specimen  test  is  assumed  representable  by  the  following 
expression 


e 

ep 


(IV-16) 


where  n  and  'P  are  temperature  dependent  material  properties,  determined,  in  the 
analysis,  from  the  pertinent  material  property  versus  temperature  curve. 


There  are  two  options  in  the  analytical  procedure  for  determining  plastic 
strains  with  the  aid  of  Equation  IV-16.  These  are  designated  as  Methods  A  and  B, 
respectively.  Method  A  involves  an  accumulative  step-by-step  process  and  is  based 
on  the  analytical  assumption  that  the  change  in  plastic  strain  induced  at  a  point  in  the 
structure  subject  to  the  stress  a',  can  be  approximated  as  Af„n  ,  shown  in 

Figure  IV-2a.  In  general,  the  plastic  strain  change  determined  in  this  manner  will  be 
an  underestimation  of  the  true  plastic  strain  change  and  is  considered  as  the  first 
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step  toward  the  computation  of  the  actual  plastic  strains.  The  plastic  strain  is  added 
to  the  previously  determined  total  plastic  strain  (  e  *  1  to  determine  new  stresses 

c  ^  for  the  same  temperature  and  load  condition.  Additional  plastic  strain  changes, 


ep 


are  computed  and  new  stresses  determined.  This  accumulative  process  is 


continued  for  either  a  specified  number  of  approximations  or  until  the  stresses  in 
successive  approximations  are  in  agreement  to  within  a  certain  number  of  significant 
digits  as  specified  by  the  analyst. 


Method  B  is  an  iterative  procedure  for  determining  the  plastic  strain  change, 
and  is  shown  in  Figure  lV-2b.  As  a  first  approximation,  the  plastic  strain  change, 

Ae  *  ,  is  determined  for  the  stress  cr*.  In  general,  the  plastic  strain  change  so 

determined  will  be  an  overestimation  of  the  true  value  in  regions  of  higher  stress  in  a 

given  structural  stress  distribution.  The  strain  change,  ,  is  then  employed 

ep 

together  with  the  previously  determined  value  (which  initially  is  zero)  to  compute  an 


which  is  used  to  obtain  the  next  approximation  to 


average  strain  change,  A  €  * 

the  stress,  or  This  iterative  process  is  continued  until  convergence  conditions 
specified  on  the  sti'esses  are  satisfied.  It  will  be  noted  that  for  Method  B  the  plastic 
changes  are  replaced  by  new  values  with  each  computational  cycle  but  for  Method  A, 
they  are  accumulated  after  each  cycle. 


The  preferred  method  to  employ  in  a  particular  plastic  analysis  will  depend 
to  a  large  extent  on  the  stress  strain  curve  and  the  details  of  the  structural  con¬ 
figuration  and  loading.  For  materials  which  exhibit  ideally  plastic  characteristics 
(i.e.  n  =  00  )  Method  A  must  be  adopted.  For  highly  redundant  or  thermally  strained 

structures,  Method  A  is  preferable,  but  for  situations  in  which  the  stresses  are  some¬ 
what  linearly  related  to  the  external  loading,  Method  B  is  recommended.  However,  it 
should  be  realized  that  cither  method  should  yield  identical  results  in  the  limit. 

The  computational  procedures  employed  to  predict  the  accumulation  of  time- 
independent  plastic  strains  will  now  be  briefly  outlined.  Consideration  is  given  to  the 
ith  time  t^  in  the  load-temperature  history  condition  under  investigation  for  a  given 

structure.  Since,  as  discussed  previously,  several  approximations  are  computed  at 
time  t.  a  particular  approximation  to  a  quantity  at  time  t.  will  be  denoted  by  the  super¬ 
script  j. 


In  Method  A,  the  steps  are  as  follows: 


(1) 


The  increment  of  plastic  strain  developed  during  the  (j-1)  th  approxi¬ 
mation  is  added  to  that  previously  accumulated  to  establish  a  total  plastic 

strain  (  e  )  pertinent  to  the  jth  approximation.  Accumulated 
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strains  associated  with  the  strain  hardening  history  are  also 
evaluated. 

J 

(2)  The  total  plastic  strains  (  €  )  arc  transformed  into  the  related 

cpT 

element  plastic  forces  {  F^1*}  and  |  }  anc'  solutions 

for  the  associated  stresses  {  a  ■*  )  and  elastic  strain  components 

(  e  '' )  are  effected  with  use  of  the  matrix  method  discussed  in  the 
e. 

1 

previous  section. 

(3)  The  signs  (tension  or  compression)  of  the  stresses  in  the  (j-l)th  and 
(j)th  approximations  are  compared  with  each  other.  If  the  sign  re¬ 
verses  the  accumulated  strain  due  to  strain  hardening  is  set  equal  to 
zero. 

(4)  Using  liquation  IV-16,  the  accumulated  strain  due  to  strain  hardening,  and 

the  elastic  strain  component,  a  series  of  operations  is  performed  which 
leads  to  the  plastic  strain  change  (Af  J  )  pertinent  to  the  jth  approxi¬ 
mation.  C*Ji 


The  result  of  step  (4)  is  employed  in  step  (1)  for  the  (j  +  l)th  approximation. 

The  process  is  continued  until  the  stresses  of  step  (2)  are  in  agreement  with  the  stresses 
of  the  previous  approximation  in  accordance  with  a  specified  criterion. 

For  Method  B,  the  computational  sequence  is  as  follows: 

(1)  The  average  increment  of  plastic  strain  (  A  e  )  determined  in  the  (j  —  1 ) th 

iteration  is  added  to  the  total  plastic  strain  as  computed  for  the  (i-l)th  time 
increment  to  define  a  new  approximation  to  the  total  plastic  strain, 


(2)  The  new  approximation  to  the  total  plastic  strain  is  then  used  in  the 
determination  of  the  (j) th  iteration  stresses,  in  the  same  manner  as  step 
(2)  of  Method  A. 

(3)  The  magnitude  of  the  plastic  strain  change  for  the  (j) th  iteration  is 
developed  with  use  of  Equation  (IV-16),  the  stresses  determined  in  (2), 

^ p  g!  i  ’  i , ; . .  i ,  iv. *  u : . .  ... ;  —  .14  4' j’ .  ...  4 1. . 


.  This  procedure  for  this  determination  differs  from  that 


of  step  (3)  of  Method  A. 
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(4)  The  (j)th  approximation  to  the  average  plastic  strain  change  is  computed 
and  the  process  of  steps  (l)-(3)  is  repeated.  The  iterative  sequence  is 
continued  until  the  specified  convergence  criteria  on  stress  is  satisfied. 

The  creep  analysis  capabilities  of  the  subject  program  are  based,  in  part,  on 
the  following  assumptions  regarding  the  mathematical  representation  of  conventional 
uniaxial  Constant  Stress-Constant  Temperature  creep  properties. 

(1)  The  primary  stage  of  the  conventional  creep  curve  is  represented  by  the 
following  relationship 


€ 

e 


O  ,  Be  ,  m 
P  (e  -l)t 


(1V-17) 


(where  e  is  the  base  for  natural  logarithms) 

whereas,  in  the  secondary  stage  one  has  the  minimum  creep  rate 
d  € 

(  e  ~  )  expressed  as 


€ 

C 


oj  (e 


B  <T 


-1) 


(IV-18) 


where  /3  ,  B,  m  and  oj  are  temperature  dependent  properties  of  the  material.  By 
equating  the  creep  rate,  as  obtained  by  differentiation  of  Equation  IV-17  with  respect 
to  time,  to  the  creep  rate  of  Equation  IV-18,  an  expression  for  the  transition  time,  t  , 
between  primary  and  secondary  ci’ecp  is  obtained  as 


(IV  -19) 


Intrinsic  in  Equation  IV-19  is  the  assumption  that  the  transition  time  is  only  a  function  of 
the  temperature  and  not  the  stress  level. 

(2)  There  exists  a  threshold  temperature,  T  ,  below  which  there  is  no 
creep  strain,  irrespective  of  the  magnitude  of  the  imposed  stresses. 

There  are  various  methods  available  for  predicting  the  accumulation  of  creep 
strains  for  the  general  condition  of  varying  stresses  and  temperatures  of  interest  here. 
The  method  adopted  for  the  present  program  is  referred  to  as  the  "strain  hardening 
rule"  which  is  fully  explained  in  Reference  9.  The  overall  computational  procedure  is 
outlined  below  for  the  ith  time  increment. 
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(1)  Computations  are  initiated  by  the  determination  of  the  average  temper¬ 
ature,  T  ,  for  the  time  increment,  from 

T.  +  T. 

T.  -  (IV-20) 

where  T.  and  T  are  the  temperatures  obtained  at  times  t  .  and  t .  respectively, 

from  die  initially  specified  temperature-time  curve.  The  average  temperature  is  next 

compared  with  the  creep  threshold  temperature,  T  .  If  T.  >  T  and  there  is  no  loss 

c  l  c 

of  creep  strain  hardening  due  to  a  stress  reversal,  there  is  a  change  in  creep  strain, 

A  e  ,  for  the  lime  increment,  A  t..  This  change  is  added  to  the  total  creep  strain 
i 

accumulated  prior  to  the  interval  to  yield  a  new  total  accumulated  creep  strain.  If,  on 
the  other  hand,  T  <  T  ,  the  change  in  creep  strain  within  the  interval  is  zero  and  the 

total  accumulated  creep  strain  from  the  previous  interval  is  carried  on  into  the  next 
interval. 


(2)  Next,  the  transition  temperature  t.  is  computed  from  Equation  I V - 1 9 

with  material  properties  obtained  for  the  average  temperature,  T..  The  type  of  creep 

(primary  or  secondary)  induced  during  the  time  increment  is  determined  by  comparing 
t.  with  an  equivalent  creep  time  parameter  t  computed  from  the  following  ex- 
1  °i-l 

pression  derived  from  Equation  IV-17  (See  Reference  9). 


'i-1 


i-1 


P 


,  ( 


i-l| 


-I 


1 

m 


(IV -21) 


where  €  c  is  the  creep  strain  accumulated  by  virtue  of  strain  hardening  (also  com- 
si-l 

puled  in  step  (1)). 

(3)  If  t  <  t.  primary  creep  prevails  and  the  magnitude  of  the  creep 
Ci-1  1 


change,  A  €  is  given  by 


A  c.  €  c. 

i  i 


(IV-22) 


i-1 


51 


ASD-TDR-63-783 

where  €  is  obtained  from  Equation  IV- 17,  written  as  follows 
c. 

1 


■<  A  t.) 


ni 


(IV-17a) 


The  creep  change  is  then  given  the  same  sign  as  the.  stress,  cr  .  If,  on  the  other  hand, 
t  >  t.  secondary  creep  is  ii.dicated  and  the  change  in  creep  strain  is  given  by 

Vl  1 


Ae.  €  At.  (IV-23) 

l  c .  l 

l 

where  the  creep  rate  e  is  obtained  front  Equation  1V-18. 

i 

(f)  Finally,  the  total  creep  strain,  e  ,  is  determined  as  thesumof  this 

°T. 

l 

change  in  creep  strain  and  the  prior  total  creep  strain  and  employed  in  the  determination 
of  the  inelastic  node  point  forces  and  a  stress  analysis  is  performed  in  the  manner  de¬ 
scribed  in  the  previous  section. 
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C.  ILLUSTRATIVE  EXAMPLES 

1.  Comparison  with  Alternate  Solutions 

The  objectives  in  presenting  the  following  illustrative  examples  are  to 

(1)  Outline  the  various  types  of  results  that  can  be  obtained  through 
exercise  of  the  respective  program  options,  and 

(2)  Demonstrate  the  level  of  accuracy  of  solutions  achieved  through  use 
of  the  program. 

To  achieve  the  second  objective  it  is  necessary  to  effect  comparisons 
with  solutions  obtained  through  alternate  techniques.  Since  these  other  techniques 
arc  limited  in  applicability  to  relatively  simple  conditions  the  following  compara¬ 
tive  discrete  element  solutions  are  necessarily  concerned  with  simple  conditions. 
The  capabilities  of  the  program  with  respect  to  irregular  geometry,  nonuniform 
load,  etc.,  are  therefore  not  demonstrated  in  these  analyses  but,  in  the  next  section, 
a  complicated  problem  for  which  only  limited  comparative  results  are  available  is 
presented. 


The  following  analyses  are  performed  in  the  present  section 

(1)  Stress  and  Deflection  Analyses 

(a)  Beam-column,  clastic 

(b)  Beam-column,  inelastic  with  initial  displacements. 

(2)  Instability  Analyses 

(a)  Elastically  restrained  column-thermal  buckling 
(Is)  Tapered  column  subjected  to  a  distributed  load. 

(c)  Column  on  four  supports. 

a.  Stress  and  Deflection  Analyses 
(1)  Beam-Column 


Figure  IV-3a  illustrates  the  conditions  for  the  analysis  of  a 
simple  beam-column  problem.  The  structure  shown  is  loaded  by  an  axial  force (P  ) 
that  is  equal  in  value  to  three-fourths  of  the  Euler  critical  load,  P  ,  i.e. 

Px=f  PE  =T  <IV-24) 


The  relationship  between  transverse  applied  loads  (in  this  case  a  concentrated  load 
(P  )  at  midspan)  and  displacement  is  linear  in  the  presence  of  a  fixed  midplane 
force.  Thus,  the  present  results  are  representable  in  nondimcnsional  form. 
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Listed  below  are  the  values  of  maximum  deflection  (deflection  at 
midspan)  and  of  maximum  bending  stress  at  midspan  as  calculated  for  idealizations 
consisting  of  various  numbers  of  elements.  Due  to  symmetry  it  was  possible  to  treat 
only  one  half  the  span.  The  parameters  employed  in  defining  these  values  were 
chosen  in  correspondence  with  the  exact  solutions,  which  can  be  written  as  follows: 

p  (sin  _v_  /V  )' (sin  J_  > 

^ - ^  ^25> 


Pz 


■yr 


y^ 


where  X 


PE 


^)(v)£ 


v  J~\  (sin  )(sin -£-^yx" 


sm 


t  yx 


(IV-26) 

(IV-27) 


and  e  is  the  distance  between  the  neutral  axis  and  the  extreme  fiber  in  compres¬ 
sion. 


Also  tabulated  below  are  the  percentages  of  error  with  respect  to 
the  exact  solution.  Figure  IV-3b  plots  the  degree  of  correspondence  of  the  calculated 
and  exact  solutions  as  a  function  of  the  number  of  elements  employed.  There  is 
extremely  close  agreement  even  for  the  crudest  idealization.  There  is  also  excellent 
agreement  between  the  deflected  shapes  computed  using  the  exact  solution  and  the 
solution  based  on  three  elements  in  the  semispan,  as  illustrated  in  Figure  lll-3c. 


No.  of 

Elements 

Per  Semispan 

Max.  Displacement 

Max.  Stress 

(r)  m 

7c  Error 

y )  (t)  (f ) 

%  Error 

1 

-0.6006 

1.72 

-3.722 

41.55 

2 

-0.6096 

0.25 

-5.437 

14.62 

3 

-0.6102 

0.15 

-5 .848 

S  .17 

4 

-0.6103 

0.13 

-6.017 

5.51 

5 

-0.6103 

0.13 

-5.106 

4.11 
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b.  Beam-Column,  Inelastic,  with  Initial  Displacements 

Sketched  below  is  a  simply  supported  axially  loaded  beam  with  initial 
displacements  (wQ)  of  the  form 

w„  0.003  sin  ZUL- 
0  20 

The  beam  is  assumed  to  be  heated  instantaneously  to  a  uniform  temperature  of  400°F. 
An  axial  load  of  3750  pounds  is  then  applied  and  sustained  at  this  value.  For  simpli¬ 
city,  the  cross-section  of  the  beam  is  assumed  to  be  composed  of  two  area  elements, 
each  equal  to  0,125  in. 2.  Since  there  is  symmetry,  it  is  only  necessary  to  consider 
half  the  beam,  idealized  as  five  elements  of  equal  length.  Numerical  values  used  for 
the  material  constants  are 


n  =  12 

cj  0.158  x  10  l> 

$ 

0.72  x  10~G  in/in 

B  = 

0.000238  in/psi 

m  =  0.72 

*  = 

0.05  x  105  psi 

E  =  9.1  x  107  psi 

a  = 

13.25  x  156  in/in/°F 

These  values  ai'e  representative  of  2024-T4  Aluminum  (See  Reference  9). 


Upon  application  of  the  axial  load,  the  beam  deflects  elastically  with 
a  total  theoretically-predicted  deflection  at  the  center  (wc)  given  by 


<"o> 


w  = 


max 


e  1  -  P/P 


Cross-Section 


P  =  3750  lb 


Deflected  Shape 


,-L 

2/3' 


20.0” 


0.125  in.-- 


2  El 

where  Pp;  =  v  -g-.  For  the  beam  under  consideration,  with  (wo)niax  =  0.005  in.. 

the  above  expression  gives  a  central  deflection  of  wc  =  0.012359  in.,  whereas  the 
program  yielded  a  value  wc  =  0.012540  in. 


The  applied  load  was  selected  small  enough  so  that  upon  its  initial 
application  only  elastic  strains  are  experienced.  With  the  passage  of  time,  however, 
creep  strains  are  sustained.  These  increase  the  deflections  and  bending  moments 
and,  therefore,  the  stresses.  This  eventually  results  in  plastic  strain. 


Results  of  the  application  of  the  computer  program  to  this  problem 
are  presented  in  Figures  IV-4  and  IV-5.  Since  the  optimum  time  increment.  At. 
was  not  known,  computations  were  performed  for  decreasing  time  increments.  Method 
B  was  employed  in  the  computation  of  time-independent  plastic  strains  and  three- 
digit  agreement  on  the  corresponding  stresses  was  specified. 


5G 


h  H 


ASD-T  DR- 03-733 


1 

1 


f 


1 

( 


[ 


Computed  total  lateral  deflections  at  the  center  of  the  beam,  presented 
as  a  function  of  time  in  Figure  IV-4a,  are  similar  to  those  obtained  in  various  exper¬ 
iments  (see,  for  example,  Reference  11).  Based  on  the  results  shown,  it  appears  that 
the  column  becomes  unstable  at  a  critical  time  of  49  hours.  As  shown  in  Figure  IV-4b, 
the  cross-sectional  area  on  the  convex  side  of  the  column  experiences  a  stress  rever¬ 
sal  at  the  critical  time.  Immediately  prior  to  the  critical  time  the  lateral  deflection 
is  finite. 

Inelastic  strains  computed  for  the  center  section  of  the  beam  are  shown 
as  a  function  of  time  in  Figure  IV-5.  As  indicated,  time-independent  plastic  strains 
are  induced  only  during  the  final  stages  of  the  deformation  process. 


b.  Instability  Analysis 

(1)  Elastically  Restrained 

The  analysis  condition  is  shown  below  in  Figure  IV-6a. 

The  beam  is  simply  supported  at  the  left  end  but  at  the  right  end  is  elastically 
restrained  by  a  torsional  spring;  transverse  and  axial  deflections  at  the  right 
end  are  prevented.  The  end  restraint  is  assigned  the  following  value. 

77  x  I01  rad/lb  in. 


e 
o 
-*— > 
o 

CJ 

QJ 

Q 


ai 


At 
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Figure  IV-4.  Total  Lateral  Deflection  at  Center  of  Beam 
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(a).  Stresses  at  Center  Section 
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(b).  Strains  at  Center  Section 


Figure  IV-5.  Axially  Loaded  Beam  -  Stress  and  Strain  History 
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The  beam  is  assumed  to  undergo  a  temperature  change  (AT)  from  the  stress  free 
state.  It  is  required  to  determine  the  value  of  AT  for  elastic  instability. 


The  element  spacing  chosen  for  analysis  is  illustrated  in  Figure 
IV-Gb;  this  spacing  was  selected  in  an  entirely  arbitrary  manner.  Based  on  these 
dimensions,  the  computations  resulted  in 


AT 

cr 


(IV -28) 


where  ci  is  coefficient  of  thermal  expansion.  By  adapting  a  solution  presented  in 
Reference  12  ,  one  can  obtain  a  comparison  solution  of 


A  T 

cr 


(I  V— 2  9) 


It  is  seen  that  the  difference  between  the  two  solutions  is  insignificant,  despite  the 
arbitrary  choice1  of  element  spacing. 

(2)  Tapered  Column  Subjected  to  Distributed  Load. 

A  problem  that  represents  a  measure  of  complexity  but  which  is 
nevertheless  solvable  in  closed  form  is  that  which  involves  the  prediction  of  the 
instability  of  a  tapered  cantilevered  column  subjected  to  a  triangular  distributed 
load  (see  Figure  IV-7a).  The  moment  of  inertia  varies  in  accordance  with  the 
expression. 

1  =  Io(1"  ~T7)  (IV-30) 


where  IQ  is  the  moment  of  inertia  at  the  support.  The  expression  for  the  varia¬ 
tion  of  "shear  load”  (q)  is 

q  =  q„  (1  -  ~  )  (IV-31) 


The  idealization  for  analysis  is  shown  in  Figure  IV-7b.  The 

distributed  loading  was  replaced  by  a  series  of  concentrated  loads  at  the  node  points 

and  for  convenience  the  reference  critical  load  (P  was  chosen  as  the  axial  load 

rel 

in  the  element  adjacent  to  the  support. 

The  significance  of  the  approximation  of  the  distributed  load  by  a 
series  of  concentrated  loads  can  be  reduced  by  developing  the  relationship  between 
qQ  and  the  axial  load  in  the  element  near  the  root.  This  relationship  is 
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7.  Buckling  of  a  Tapered  Column  Subjected 
a  Triangular  Distributed  Axial  Load 
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( 


ii 
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) 


cr 


1 


(3N-4) 

3(N-1)2 


where  N  is  the  number  of  node  points  in  the  idealization.  Thus,  the  evaluated 
critical  reference  loads  are  divided  in  the  indicated  manner  to  obtain  the  calcualted 
critical  value  of  distributed  load.  For  purposes  of  comparison  with  the  exact  solu¬ 
tion,  the  results  of  analyses  for  various  numbers  of  node  points  are  shown  below. 
For  purposes  of  comparison  with  the  exact  solution  as  given  in  Reference  12  , 
pg  139,  the  results  are  presented  in  the  form  of  the  ratio  between  fio^/2)  and  the 

2 

Euler  critical  load  (P ...  =  t t  ~EI0  )  of  a  simply  supported  beam  of  the  same  length 

L2 

and  moment  of  inertia  I0.  The  "exact"  value  of  this  ratio  is  1.317.  Figure  IV-7c 
gives  a  graphical  representation  of  the  results.  Note  that  in  contrast  with  uniform 
load  and  geometry  conditions,  the  present  case  requires  10  node  points  to  reduce 
the  error  below  engineering  significance. 


Number  of 
Node  Points 

q0L/  (calculated) 

— /PK 

’/<_  Error* 

3 

0.912 

<30.7 

-1 

1.117 

<  15.2 

5 

1.201 

-t  9.7 

11 

1.299 

+  1.-11 

15 

1 .308 

+  0.7 

*%  Error  = 


1.317 


Wp 

2  '  E 


1.317 


(3)  Column  on  Four  Supports 

As  an  example  of  a  more  complicated  instability  problem,  Die 
case  of  the  uniform  section  column,  supported  at  four  points,  has  been  solved.  The 
conditions  of  analysis  arc  illustrated  in  Figure  IV-8.  The  outer  ends  are  simply 
supported  while  the  two  inner  supports  are  elastically  restrained.  Solutions  for  this 
problem  were  developed  by  Klemperer  and  Gibbons  (Reference  13)  and  described  by 
Timoshenko  in  Reference  10,  pg.  107. 

In  performing  the  present  analyses,  two  gridworks  were  employed: 
a  subdivision  of  each  segment  between  supports  into  three  elements,  and  also  a  sub¬ 
division  into  four  elements.  The  rigidity  of  the  supports  is  defined  by  the  restraint 
coefficient,  (w.  However,  for  purposes  of  data  representation  it  is  preferable  to 
embody  this  consideration  into  a  "support  rigidity  parameter",  defined  as  (  L  , 
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where  P_  is  tlie  Kuler  critical  load  for  a  single  span,  P  - —  .  The  results, 

E  E  1/ 

expressed  as  the  ratio  between  the  critical  value  of  the  applied  load,  P(,r,  and  P  ai'e 

tabulated  below  for  a  wide  range  of  support  rigidity  parameters  and  are  plotted  in 

Figure  IV-8b. 


Figure  IV-Sb  also  shows  the  results  presented  in  Reference  10  , 

Figure  72.  There  is  extremely  close  agreement  between  the  present  results  and 

those  given  in  Reference  10.  It  is  noteworthy  that  three  distinct  relationships 

between  the  parameters,  P  / P ,./  ,  and  C  L/P.,  ,  governed  by  three  different 

e  i'  h  w  E 

mode  shapes,  arc  involved.  These  mode  shapes,  as  determined  by  the  output  of  the 
present  analyses,  are  shown  in  Figure  lV-8c;  their  regions  of  occitrance,  as  defined 
in  Reference  10,  are  shown  in  Figure  lV-8c.  It  is  noteworthy  that,  for  values  of 
support  rigidity  parameters  that  lie  within  the  center  region  and  near  the  curve 
intersections  it  was  difficult  to  obtain  convergence  to  the  lowest  root.  Support 
rigidity  parameter  values  in  excess  of  80  required  over  150  iterations  for  conver¬ 
gence. 
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2,  Analysis  of  Practical  Complex  Conditions 


As  an  example  of  the  analysis  of  a  complex  indeterminate  one-dimensional 
structure  for  thermal  stress  and  applied  load,  the  conditions  illustrated  in  Figure 
IV-9a  are  examined.  The  conditions  shown  were  analyzed  for  support  reactions  in 
Reference  14.  The  present  analysis  considers  not  only  the  equilibrium  solution  for 
stress  and  displacement  but  also  solves  for  one  form  of  elastic  instability. 

The  discrete  element  idealization  employed  appears  in  Figure  IV-9  b  and 
c.  The  entire  length  of  the  structure  is  divided  into  nine  elements  each  of  which  is 
assigned  10  differential  areas.  In  addition,  there  is  a  flexible  axial  support  at  the 
right  end.  The  material  properties  are  specified  as:  K  =  107  psi,  and  a  10-5  in./ 
in.°F,  while  the  temperature  distribution  is 


T  =  2x  -t  (0.4x2  -  2X)£  0  <  x  <  10  (IV- 33) 

T  =  20  <  20f  10  <  x  <  20  (IV -34) 

(It  is  assumed  that  the  temperature  in  the  stress-free  state,  TQ,  is  0°F).  Thus,  the 
temperatures  vary  both  in  the  axial  and  deptlnvise  directions.  In  the  left  span,  the 
temperatures  on  each  differential  area  of  an  element  take  on  the  value  at  the  center 
of  the  element. 

—  Results  arc  presented  in  Figure  IV-10.  The  support  reactions  are  in  close 
correspondence  with  the  results  of  Reference  14  (see  Figure  IV-lOa).  The  differences 
between  the  two  sets  of  values  arc  due  to  the  fact  that  the  Reference  14  results  ai-e 
obtained  from  an  essentially  "closed  form"  solution  while  the  present  values  include 
discretization  errors.  The  variation  of  extreme  fiber  stresses  associated  with  these 
reaction  forces  is  also  shown  in  Figure  IV-lOa. 

Figure  IV-lOb  shows  the  displacement  patterns  due  to  temperature  alone, 
temperature  and  applied  load  combined,  and  buckling. 
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Figure  IV-9.  Illustrative  Complex  Indeterminate 
Beam  Problem 
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CHAPTER  V 

PLATE  ANALYSIS  PROGRAM 


A.  SCOPE 

The  objective  of  the  plate  analysis  program  is  to  determine  the  structural 
behavior  (i.e.,  stresses,  displacements,  and  buckling  stress)  of  flat  plates  of  arbi¬ 
trary  planlorm  and  thickness  variation  subjected  to  nonuniformly  distributed  applied 
loads  and  temperatures.  The  plate  material  can  be  orthotropic  but  the  principal 
directions  of  ortholropy  and  corresponding  material  properties  must  be  the  same  at 
all  points  throughout  the  plate.  Stiffening  members  can  be  accommodated;  the  pro¬ 
gram  can  in  fact  be  used  for  truss  or  beam-gridwork  analysis.  The  temperature 
dependence  of  material  properties  as  well  as  inelastic  behavior  can  be  taken  into 
account. 

Three  types  of  discrete  elements  are  accommodated  in  the  program: 

(1)  The  beam  segment  (Fig.  Ill- 2). 

(2)  The  triangular  plate  (Fig.  Ill-  (i). 

(3)  The  quadrilateral  plate  (Fig.  111-7  ). 

The  behavior  of  each  of  these  elements  is  segregated  into  ''inplane"  (or  axial)  be¬ 
havior  and  "out-of-plane”  (or  flexural)  behavior.  Relationships  were  formulated 
in  Chapter  III  for  the  beam  segment.  The  bases  for  the  derivation  of  the  properites 
of  plate  elements  are  discussed  in  Chapter  III;  details  are  given  in  Reference  8. 

Figure  V-l  shows  a  typical  plate  structure  whose  behavior  is  to  be  predicted  by 
the  use  of  the  program.  Assume  further  that  it  possesses  orthotropic  material 
properties  with  the  principal  directions  of  ortholropy  being  parallel  to  the  x  and  y 
axes.  The  reference  axes  for  analysis  must  correspond  to  the  principal  directions 
of  ortholropy. 

Once  the  system  axes  have  been  established,  an  arrangement  of  discrete  ele¬ 
ments  in  idealization  of  the  actual  structure  must  be  decided  upon.  Program  capa¬ 
city  limits  the  scope  of  the  idealization  to  a  maximum  of  80  node  points. 

If  the  thickness  of  the  actual  structure  varies  continuously,  this  variation  must 
be  represented  in  a  stepped  manner.  With  regard  to  the  choice  between  the  triangular 
and  quadrilateral  elements,  it  is  believed  that  the  quadrilateral  is  sufficiently  versa¬ 
tile  for  most  applications  and,  in  addition,  yields  a  solution  for  stress  that  is  more 
easily  interpreted  than  the  solution  for  stress  in  a  triangular  plate.  The  triangular 
plate  is  indispensable,  however,  under  certain  geometric  conditions. 
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(d)  Section  A-A  —  Idealization 


Figure  V-l.  Typical  Plate  Analysis  Problem 
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Applied  direct  loads  can  exist  in  either  the  x,  y,  or  z  directions  at  each  node 
point,  and  applied  moments  about  the  x  and  y  axes  can  be  present  at  each  node. 
Normally,  loads  will  exist  only  at  the  node  points  along  the  edges  of  the  structure  and 
they  will  likely  be  initially  defined  as  distributed  stresses.  It  is  the  responsibility  of 
the  analyst  to  transform  the  distributed  applied  stresses  into  equivalent  concentrated 
forces  at  the  affected  node  points.  Subject  to  previously  cited  limitations,  as  many  as 
10  independent  applied  load  conditions  can  be  solved  for  in  a  single  computational 
sequence . 


The  temperature  distribution  is  assumed  to  be  known  at  the  start  of  the  analy¬ 
sis.  It  must  be  emphasized  that  only  one  temperature  condition  per  computational 
sequence  is  permissible  since  the  stiffness  properties  of  the  system  are  dependent 
upon  the  material  mechanical  properties  which,  in  turn,  are  allowed  to  be  dependent 
upon  temperature.  In  general,  a  change  in  temperature  conditions  requires  a  recal¬ 
culation  of  the  equations  governing  the  elastic  behavior  of  the  structure. 


The  temperature  state  of  the  idealization  for  inplane  analysis  is  initially 
defined  by  the  temperatures  of  the  node  points.  During  the  computational  process 
each  element  is  assigned  a  uniform  temperature  value  which  is  the  simple  average  of 
the  temperatures  at  its  corner  points. 


In  the  plate  idealization  for  out-of-plane  analysis,  the  input  must  define  the 
'thermal  - -  '~*a 


moments"  (M  ,  M  )  at  the  node  points.  These  are 
lL  ■' 

Ey(l  +  *  xy) 


(1  ~  u  u  ) 
'  xy  '  yx' 


E  (1  +  u  ) 
x  'yx  > 


<1- 


11  xy  'V 


a  T  i  d  £ 


oTl  dC 


(V-l) 


Temp. 


where  {  is  a  coordinate  measured  normal  to  the  middle  surface.  For  the  particular 
case  of  an  isotropic  plate  of  thickness  t  with  a  linear  temperature  gradient  through 
the  thickness  and  constant  material  properties 


_  a 
M 


E  a  (Ti-T2)h2 

12(1-  /i.  ) 


(V  — 2) 


where  T^  and  T^  are  the  outer  and  inner  surface  temperatures,  respectively.  During 
the  computational  process,  each  plate  element  is  assigned  a  thermal  moment  which  is 
the  simple  average  of  the  thermal  moments  at  its  corner  points. 
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For  each  flexural  clement,  the  thermal  moment  at  each  end  (node  point)  must 
be  included  with  the  input.  The  formula  for  computing  this  value  is: 

|  i  Temp. 

rA  i  l  'V7?T 

M.  «  /  EOT;  cl  A  |  -  r  H - /// 


|  Elastic  Axis 

or,  if  the  material  properties  arc  regarded  as  constant  for  the  purpose  of  this  deter¬ 
mination  and  the  cross-section  is  a  solid  rectangle,  then 


E  q  (Ti-T2)  b  d~ 


Temp. 


The  material  mechanical  properties  E  ,  E  ,  u  ,  u  ,  G  ,  and  the 

x  y’  xy’  yx’  xy 

coefficient  of  thermal  expansion  (  a)  are  each  permitted  to  sustain  an  independent 
variation  with  temperature,  as  represented  by  as  many  as  five  points  on  the  material 
property  versus  temperature  curve.  This  would  appear  to  be  a  sufficient  number  of 
points  for  the  common  sti'uctural  materials.  If  the  true  relationship  is  extremely 
complicated  across  the  full  range  of  temperatures  it  is  likely  that  the  problem 
involves  only  a  restricted  portion  of  this  range,  which  in  itself  can  be  represented  by 
five  points.  When  evaluating  the  material  properties  for  a  given  element  the  compu¬ 
ter  first  selects  the  element  temperature  and  then  establishes  the  desired  properties 
by  linear  interpolation  of  the  material  property  versus  temperature  data. 

Capabilities  for  inelastic  analysis  included  in  this  program  relate  only  to  time 
independent  behavior.  The  basic  terms  for  inelastic  analysis  are  included  in  the 
next  section,  which  is  a  review  of  the  computational  process  for  elastic  analysis. 

As  in  the  case  of  the  beam  program  (Chapter  IV)  inelastic  analyses  arc  performed 
as  a  series  of  elastic  analyses.  The  incorporated  rules  lor  material  inelastic 
behavior  as  well  as  the  procedures  which  govern  the  series  of  analyses  for  such 
behavior  arc  examined  in  a  later  section. 

B.  THEORETICAL  BASIS 

1.  Elastic  Analysis 


The  following  discussion  of  the  analysis  procedure  begins  with  a  treat¬ 
ment  of  the  inplanc  analysis  sequence  since  these  operations  arc  always  performed 
first  by  the  program.  The1  relationships  between  the  inplane  corner  point 
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displacements  (u,  v)  of  each  element  and  the  corresponding  corner  point  forces  are 
expressed  in  the  form: 

&  w  « -{?}-{?} 

where  [k  1  is  the  element  "inplane"  stiffness  matrix,  ■[  F  ,  F  ]■  are 
L  xyJ  f  p  pi  <■  *  y  J 


"thermal  forces"  at  the  node  points  and  j  F^’F^Jarc  the  node  point  plastic  forces 

(assumed  known  in  this  discussion).  The  algebraic  form  of  these  relationships  is 
detailed  in  Chapter  III. 

Based  on  the  pertinent  input  data  —  geometry,  loads,  temperatures,  and 
material  properties  —  the  computer  first  evaluates  the  inplane  element  relationships 
and  constructs  with  them  the  system  of  equations  that  represent  the  analytical  model 
of  the  complete  structure  for  inplane  analysis.  These  assembled  equations  are  of 
the  form 


KvI  C} 


where 


F  1 

x  r  i 

“  are  applied  loads  at  the  node  points,  K  is  the  "master  inplane 

y/ .  ..  ,fpQ1  .  fppl  ..  J..xy;  . 


x  >are  the  "net"  thermal  and  plastic  forces 


stiffness  matrix",  and 


at  the  node  points.  Although  Equation  V-6  indicates  a  single  column  of  applied  loads, 
as  many  as  10  different  applied  load  conditions  can  be  treated  in  an  inplane  analysis 
cycle  for  a  given  temperature  distribution  and  in  the  absence  of  inelastic  behavior. 
Only  one  temperature  condition  can  be  treated  since  the  stiffness  matrix  j^K^J  is 

a  function  of  the  material  properties  and  these  in  turn  are  a  function  of  temperature. 


For  each  designated  support  condition  a  column  and  the  corresponding  row  are 
eliminated  from  the  [l<  ]  matrix  and  the  affected  rows  are  removed  from  the 

r  Pv  vT  , r  „  en  p  i 


column  matrices 


,  and  J  Px‘ 


Utilizing  the  sub¬ 


script  "R"  to  designate  the  thusly  reduced  matrices,  one  obtains 


=  r  k  { u} 

l  xyj  „  lvJ„ 


(V-Ga) 
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The  basic  unknowns  in  Equation  (V-Ga)  are  the  displacements 
solved  for  as  follows: 


and  these  are 


Once  the  solution  for  the  displacements  has  been  achieved,  a  "listing"  is  made 
of  the  displacements  associated  with  the  respective  elements,  so  that  they  can  be 
applied  to  the  determination  of  the  element  stresses.  Relationships  between  the 
stresses  in  the  element  and  the  element  node  point  displacements  are  also  stored  in 
the  program  in  the  following  form: 


(V-8) 


jer^,  T  ,  0-y  |  is  a  listing  of  the  stresses  at  the  corner  points  of  the  element, 

[S  ]  is  the  "element  stress  matrix",  and  j  <7^7  0,  '  |  represents  the 

stresses  developed  on  the  basis  of  inducing  the  thermal  expansion  elastically. 
Equations  V-8  are  discussed  in  Chapter  III. 


The  procedure  employed  for  out-of-plane  analysis  is  quite  similar  to  that 
described  above  except  that  in  the  presence  of  instability  effects  an  additional  or 
"incremental"  stiffness  results  from  the  presence  of  midplane  forces.  The  element 
equation  is  now  written  as 


{ 


M 

x 


} 


> 


(V-  9) 


where  j^nj  is  the  incremental  stiffness,  the  terms  of  which  are  functions  of  the 
midplane  stresses  determined  through  Equation  (V-9) .  The  detailed  form  of 
Equation  (V- 9)  is  examined  in  Chapter  III.  Again,  the  individual  element  stiffness 
matrices  are  evaluated  and  summed  to  form  the  stiffness  matrix  for  the  complete 
structure  and  the  support  conditions  are  applied,  resulting  in 
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(V-10) 


It  "  R  R 

If  the  solution  for  the  equilibrium  displacements  and  stresses  resulting  from 
applied  loads  and  temperature  change  is  sought,  one  obtains: 
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^(V-ll) 


and,  from  the  appropriate  relationships 
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(V-12) 


It  is  to  be  noted  that  moments  and  shears  per  linear  inch,  rather  than  stresses,  are 
employed  to  define  the  internal  distribution  of  load.  (Such  moments  are  "primed" 
rather  than  "barred"  to  indicate  a  distinction  with  the  concentrated  corner  point 
moments). 

If  critical  stresses  are  to  bo  determined  (i.e.,  if  an  elastic  instability 


analysis  is  to  be  performed)  j  M^,  M^,  |  and  |  M ”,  M°,  p7°} 


R 


are  set  equal  to  zei-o,  the  matrix  N  J  is  multiplied  by  the  scalar  X  ,  and  Equation 
V-10  is  rearranged  as  follows 


X 


-1 


r  k  l 

r  n] 

l  zj  R 

L  J 

It 


(V-9) 


The  scalar  X  (which  is  the  eigenvalue  to  be  determined)  as  well  as  the  relative 
magnitude  of  the  displacements  j  0,,  0  ,  w  j.  (the  eigenvector,  which  is  the 

buckled  shape)  are  then  determined  through  matrix  iteration  as  described  in  Chapter 
II,  Note  that  X  represents  a  value  by  which  all  applied  midplane  loads  are 
multiplied  to  achieve  instability.  Thus,  if  X  =  1.0,  instability  is  reached  and  if 
X  >1.0  the  given  midplane  loads  arc  not  sufficient  to  produce  instability. 

The  matrix  £  N  ]  is  dependent  upon  the  values  of  stress  computed  in  the 
inplane  analysis  portion  of  the  computational  cycle.  If  many  load  conditions  are 
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cited,  the  program  will  automatically  utilize  the  inplane  stresses  due  to  the  first 
inplane  load  condition  in  the  formulation  of  the  [n  j  matrix. 

An  additional  capability  of  the  program  is  its  ability  to  treat  "oblique" 
support  conditions.  It  sometimes  occurs  that  a  structure  is  constrained  to  displace 
in  the  direction  of  axes  other  than  those  employed  in  the  definition  of  the  behavior 
of  the  structure  as  a  whole.  (See  sketch).  Such  conditions  can 


be  accommodated  simply  by  specifying  for  each  affected  point,  the  direction  of  the 
special  x'  axis  by  means  of  node  points.  The  program  will  then  transform  the 
relationships  at  all  these  points  into  equations  referenced  to  the  special  axes,  and 
boundary  (support)  conditions  can  be  defined  with  respect  to  the  latter. 

2.  Inelastic  Analysis  Procedure 

An  extension  to  the  above-described  plate  analysis  program  permits 
analyses  for  time-independent  plastic  behavior  under  conditions  of  varying  stress 
and  temperature.  This  extension  is  limited  to  the  use  of  the  triangular  plate  element 
and  excludes  orthotropic  behavior,  i.c.,  only  isotropic  behavior  can  be  treated. 

Consider  the  typical  triangular  plate  element  sketched  below.  Since  the 
stress,  and  therefore  the  plastic  strain,  varies  across  the  thickness  of  the  plate,  it 
is  necessary  for  purposes  of  plastic  analysis  to  divide  the  thickness  into  laminae, 
or  sub-elements,  of  thickness  A  t.  The  location  of  a  given  sub-element  is  specified 
by  the  £  value  to  the  center  of  the  element,  as  shown.  As  many  as  10  sub-elements 
can  be  employed  in  a  given  analysis. 
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To  define  the  temperature  history  of  the  structure  it  is  necessary  to 
specify  the  initial  temperature  of  the  entire  plate  and  a  temperature  history  for  each 
lamina  at  every  node  point. 

The  effective  temperature  for  each  sub-element  is  then  computed  as  the  arithmetic 
mean  of  its  three  corner  temperatures.  These  average  temperatures  are  used  to 
determine  the  pertinent  material  properties  and  the  temperature  differences  A  T. 


At  each  specified  time  in  the  load-temperature  history,  equivalent  thermal 
stress  resultants  for  inplane  and  flexural  behavior,  as  required  for  construction  of 
the  comer  thermal  forces  and  moments,  are  computed  using  the  following  equations 


Na  =  - - - 

(l  -H-  ) 

X  E  a  A  T  A  h 

(V-13) 

M  °  -  1 

X  Ea  AT  Ah 

(V-14) 

(1  -  M  ) 

Note  that  Equation  (V-14)  is 

a  modification  of  Equation  (V-l). 

Note  also  that  in  the 

inelastic  analysis  routine  the  thermal  moments  are  evaluated  by  the  computer, 
based  on  input  data,  while  for  elastic  analyses  these  moments  must  be  hand  com¬ 
puted  and  entered  with  the  input. 


In  addition,  the  inplane  and  flexural  stiffnesses  required  in  the  deter¬ 
mination  of  the  elastic  neutral  axis  £  and  in  the  construction  of  the  stiffness 
matrices  are  given  by 


Eh  =  £  E  A  h 

El  =  <£  ')2  A  h 


where 


£ 


£  -  £ 


and 


r=  £e  i  A  h 

^  Eh 


(V-15) 

(V-16) 


The  various  inelastic  column  force  matrices  present  in  the  general  force 
displacement  relationships  are  composed  of  the  following  inelastic  stress  resultants. 

Nx  =  “J-2-IE(exP  +  />  Ah  MP=  _L^XE(exP+  ^%P)£  ‘  Ah 

(1-/0  XT  yT  X  (1  -M2)  XT  ^T  ' 


N"  = - ^~XE(€  ,P  +  H-*  ,P  )  Ah 


yT  *T 


XE(e/  +  /*«  F>  £  '  Ah 

(1-0  yT  T 


M  =  „ 

y  ..  2 


N 


xy  2(1+0 


XttXe  y  p  a  h 


M  P  =  XE  /  P  £  '  Ah 

xy  2(1  +  H-  )  ^  xyT  * 

(V- 17) 
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The  symbols  e  ^  ,  e  and  Y  ^  represent  the  total  plastic  strains  for 
-nt  yT  xy-p 

the  biaxial  stress  case  and  are  computed  with  the  aid  of  plastic  analysis  Method  A 
discussed  in  Chapter  IV  for  the  uniaxial  stress  case.  By  consideration  of  the 
uniaxial  stresses  and  strains  employed  in  Method  A  as  effective  stresses  and 
strains,  and  adopting  the  biaxial  inelastic  stress-strain  relationships  consistent  with 
the  incremental  theory  of  plastic  flow,  plastic  strain  increments  /\  e  P  , 

XT 

A  6  ^  ,  /  P  a  re  computed  and  accumulated  for  a  given  load- temperature 

yT  xy 

history.  Details  of  this  portion  of  the  procedure  are  delineated  in  Reference  9. 

The  extension  of  the  plate  pi'ogram  requires  the  detailed  computation  of 
the  stresses  for  each  area  from  the  following  expression 
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( V —18) 


The  general  computational  procedure  parallels  that  of  the  one-dimensional 
program.  First,  the  load-temperature  history  is  divided  into  time  intervals  and 
each  of  these  is  in  turn  subdivided  into  time  increments.  Consider,  for  example,  the 
computation  of  the  stresses  at  time  t.  .  At  this  time  the  elements  of  the  plate  may 
have  experienced  prior  biaxial  plastic  straining.  Node  point  deformations  are  com¬ 
puted  for  the  imposed  loading  and  net  plastic  forces  and  are  then  employed  in  the 
computation  of  the  biaxial  stresses  from  the  element  stress  equations.  These 
stresses  are  used  to  evaluate  effective  stresses  which  are  used  to  compute  the 
change  in  effective  strains  requii'ed  subsequently  in  the  determination  of  the  change 
in  biaxial  plastic  strains.  The  latter  are  added  to  the  prior  total  plastic  strains  to 
establish  new  values  for  the  total  plastic  strains.  Stresses  arc  then  computed  and 
the  process  continued  until  successive  stresses  ai'C  reproduced  in  satisfaction  of  a 
convergence  criteria. 
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C.  ILLUSTRATIVE  EXAMPLES 

1.  Comparison  with  Alternate  Solutions 

As  in  the  case  of  the  one-dimensional  structure  program,  the  objectives 
in  presenting  the  following  illustrative  examples  are  to  provide  an  outline  of  the 
results  to  be  obtained  through  exercising  the  various  program  options  and  to  demon¬ 
strate  the  level  of  accuracy  that  can  be  achieved  with  use  of  the  program.  To 
accomplish  the  second  objective  it  is  again  necessary  to  compare  results  with  solu¬ 
tions  for  relatively  simple  conditions.  For  the  present  program,  the  comparisons 
are  further  limited  in  versatility  by  the  relatively  few  available  alternate  solutions 
for  orthotropic  plate  problems.  Although  the  theoretical  formulation  of  the  governing 
differential  equations  for  orthotropic  behavior  are  well  established,  there  are  few 
numerical  solutions  to  such  problems.  No  significant  solutions  for  thermal  stress 
conditions  in  orthotropic  plates  have  been  found  by  the  authors.  Hence,  in  what 
follows,  the  thermal  stress  problems  are  concerned  only  with  isotropic  plates. 

The  following  types  of  analyses  are  performed  in  this  section; 

(a)  Stress  and  Deflection  Analyses 

(1)  Isotropic  Rectangular  Plate  -  Thermal  stress  analysis 

(2)  Orthotropic  Rectangular  Plate  -  Plane  sti'ess  analysis 

(3)  Isotropic  Triangular  Plate  -  Thermal  sti'ess  and  flexure 

(b)  Instability  Analyses 

(1)  Isotropic  Rectangular  Plate  -  Thermal  buckling 

(2)  Orthotropic  Rectangular  Plate  -  Flexure  and  Instability 

Analyses 

The  analysis  of  a  problem  involving  geometric  irregularities  and  more  practical 
circumstances  than  the  above  problems  will  be  described  in  the  next  section. 


a.  Stress  and  Deflection  Analyses 

(1)  Isotropic  Rectangular  Plate  -  Thermal  Stress  Analysis 

The  rectangular  plate  shown  in  Figure  V-2a  has  a  length  that  is 
twice  its  width  and  has  a  thickness  variation  in  the  width  direction  given  bv 
h  =  ["1-0.9 

2 

T  =  [(“)  ~  "3  ]  Tq,  producing  midplanc  clastic  thermal  stresses.  A  solution  for 

these  thermal  stresses  was  advanced  by  Mendelson  and  Ilirschbcrg  in  Reference  15  . 


y  .  ^  J  hQ  .  The  temperature  varies  in  the  width  direction, 
c  ' 


79 


ASD-TDR-63-783 


Temperature  Distribution: 


Figure  V-2.  Isotropic  Rectangular  Plate  for  Thermal  Stress  Analysis 
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b.  Idealization  of  Heated  Plate  -  One  Quadrant 


Figure  V-2.  (Concl'd)  Isotropic  Rectangular  Plate  for  Thermal  Stress  Analysis 
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From  symmetry  it  is  only  necessary  to  consider  a  quadrant  of 
the  plate.  The  idealization  for  analysis  consists  in  dividing  tine  plate  quadrant  into 
54  quadrilateral  plate  elements  as  illustrated  in  Figure  V-2b.  Judgement  dictates 
the  placement  of  small  elements  at  the  tip  and  side  where  stresses  change  rapidly  in 
consequence  of  edge  effects  and  the  temperature  gradient  and  thickness  changes 
are  most  severe. 


Selected  results  from  the  present  analysis  and  from  Reference 
15  are  plotted  in  Figure  V-3.  All  results  are  expressed  in  nondimensional  form. 

The  distribution  of  longitudinal  thermal  stress  (  <r is  given  on  a  cross-section 

located  0.25c  from  the  end  (x  =  1.75c),  where  this  component  of  stress  is  signifi¬ 
cantly  influenced  by  "end  effects".  The  chordwise  stress  (  a  ),  which  is  entirely 

due  to  end  effects,  is  plotted  for  the  tip  chord  (x  =  2.00c)  while  the  chordwise 

variation  of  shear  stress  (  T  )  is  shown  for  x  =  1.60c.  The  dotted  lines  represent 

xy 

the  computer  program  solution;  the  solution  from  Reference  15is  shown  by  solid 
lines.  As  seen  from  this  figure,  there  is  close  agreement  between  the  two  solutions. 

(2)  Orthotropic  Rectangular  Plate  -  Plane  Stress  Analysis 

For  plane  stress,  one  of  the  most  interesting  orthotropic  plate 
problems  to  have  been  solved  is  illustrated  in  Figure  V-4a.  An  orthotropic  plate, 

E  ,  =  0.17  and  G_  =  0.134,  of  aspect  ratio  4.0,isloaded  by  equal  and  opposite 

y/  X  Ey 

concentrated  forces  P  in  the  manner  indicated.  The  plate  thickness  is  h. 

The  idealization  for  discrete  element  analysis  consists  of  the  36 
rectangular  plate  elements  shown  in  Figure  V-4b  (due  to  symmetry,  only  a  quadrant 
of  the  plate  need  be  treated).  Again,  smaller  elements  are  utilized  in  the  vicinity  of 
the  concentrated  load.  The  discrete  clement  analysis  results  are  presented  in 
Figure  V-4c,  where  they  are  compared  with  the  results  derived  by  Conway 
(Reference  16  )  by  use  of  a  classical  type  of  approach.  It  is  seen  that  there  is  an 
extremely  close  agreement  between  the  two  solutions. 

(3)  Isotropic  Triangular  Plate  -  Thermal  Stress  and  Flexure 

To  illustrate  the  accuracy  of  the  program  in  the  analysis  of  plates 
subjected  to  temperature  gradients  across  the  thickness,  the  equilateral  triangular 
plate  shown  in  Figure  V-5a  has  been  analyzed.  The  linear  450 °F  temperature 
gradient  through  the  thickness  is  assumed  to  be  constant  throughout  the  plate.  An 
analytical  solution  for  this  type  of  problem  was  obtained  by  Maulbetsch  (Reference 
17  ). 


Since  the  plate  is  symmetric  about  the  x-axis  it  was  necessary  to 
consider  only  one  half  of  the  plate.  The  idealization  was  achieved  by  means  of  21 
quadrilateral  plate  elements  and  7  triangular  elements  as  shown  in  Figure  V-5b. 
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(a)  Actual  Plate 


(b)  Analytical  Idealization  —  One  Quadrant 


-0.13GP  0.241P  0.259P  0.294P  0.370P 


— Stress  (  ^y) 


0.267P  2.340P  1.375P  0.740P  0.615P 
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(c)  Comparison  of  Solutions  for  Longitudinal  Stress 

-  Discrete  Element  Idealization  Results 

- Alternate  Analytical  Approach 


Figure  V-4.  Inplane  Stress  Analysis  —  Rectangular  Orthotropic  Plate 
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(a)  Equilateral  Triangular  Plate 


Figure  V-5. 


Isotropic  Triangular  Plate  for  Thermal  and  Flexural  Analysis 
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To  satisfy  the  support  conditions  along  the  swept  edge,  it  was  convenient  to  apply 
"oblique"  coordinates  along  this  edge. 

The  results  of  the  discrete  element  analysis  are  compared  to  the 
results  of  Reference  17  in  Figure  V-6.  The  deflection  of  the  plate  along  the  x-axis 
is  plotted  in  Figure  V-6a.  An  almost  exact  agreement  exists  between  the  deflec¬ 
tions  obtained  from  the  program  and  the  analytical  results  of  Reference  17. 

The  thermal  moments  Mx  and  My  (moments  about  the  x-axis 
and  y-axis  respectively)  are  independent  of  the  y-coordinate.  The  distribution  of 
these  moments  in  the  x-direction  are  shown  in  Figure  V-6b.  The  symbolized  points 
for  the  program  results  correspond  to  the  elements  adjacent  to  the  x-axis.  With 
the  exception  of  two  elements,  the  results  are  in  very  good  agreement. 

Figure  V-6c  compares  the  thermal  twisting  moments  MXy  by 
plotting  their  distribution  in  the  y-direction  (the  twisting  moments  are  independent 
of  the  x-coordinate) .  The  discrete  element  results  given  in  Figure  V-Gc  correspond 
to  the  elements  adjacent  to  the  y-axis.  The  results  are  again  in  very  good  agree¬ 
ment  except  for  the  triangular  and  adjacent  rectangular  element  which  are  in  fair 
agreement. 

b.  Instability  Analyses 

(1)  Isotropic  Rectangular  Plate  -  Thermal  Buckling 

As  an  illustration  of  buckling  caused  by  temperature  change, 
the  problem  of  a  uniformly  heated  rectangular  plate,  simply  supported  along  two 
parallel  edges  and  restrained  against  rotation  and  normal  expansion  along  the 
other  two  edges,  is  analyzed.  The  conditions  of  analysis  are  as  shown  in  Figure 
V-7  where  4  rectangular  elements  have  been  employed  in  the  idealization  of  a  quad¬ 
rant  of  the  plate. 


The  theoretical  buckling  stress  in  the  x-direction  for  this  plate 
is  determined  by  the  expression 
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Using  the  value  of  the  buckling  coefficient  k  as  6.65  (from  Reference  10),  the  criti- 

.  2  C 

cal  stress  is  computed  to  be  4050  lb/in.  . 

2  A  temperature  rise  of  15.9°F,  corresponding  to  a  thermal  stress 
of  2000  Ib/in.  ,  was  imposed  on  the  plate  and  an  inplane  and  instability  analysis 


8G 


Twisting  Moments  Thermal  Moments  ~  Deflection. 

~  103  in.-lb/in.  103  in.-lb/in.  w,~  inches 


ASD-TDR-63-783 


Figure  V-G.  Comparison  of  Results  -  Isotropic  Triangular  Plate 
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(a)  Plate  Planform  and  Support  Conditions 


(b)  Idealization  of  Plate  Quadrant 

Figure  V-7.  Isotropic  Plate  for  Thermal  Buckling  Analysis 
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performed  by  the  subject  computational  program.  The  results  of  the  inplanc  analy¬ 
sis  indicated  a  compressive  stress  of  2003  lb/in.2  in  the  x-direetion  and  the 
instability  analysis  showed  a  critical  inplanc  load  factor  (  — £-)  of  1.9015.  The  pro¬ 
gram,  therefore,  predicted  a  buckling  stress  of  1.9015(2003)  3860  lb/in.2  which 

agrees  within  5  percent  of  the  theoretical  value  of  4050  lb/in.2.  The  accuracy  of  the 
predicted  buckling  stress  may  be  improved  by  utilizing  a  larger  number  of  elements 
in  the  idealization. 

(3)  Orthotropic  Rectangular  Plate  -  Flexure  and  Instability  Analysis 

Formulas  for  the  deflection  of  rectangular  orthotropic  plates 
subjected  to  a  concentrated  lateral  center  load  and  uniform  lateral  pressure  have 
been  published  by  liearmon  (Reference  18),  and  formulas  for  the  instability  of 
orthotropic  plates  under  various  support  conditions  are  presented  in  Reference  19. 
The  following  is  a  comparison  of  certain  of  these  formulas  with  results  obtained 
by  discrete  element  analysis  for  the  particular  conditions  illustrated  in  Figure  V-8a. 
Two  cases  are  examined: 

(1)  The  deflection  of  the  center  of  the  plate  under  a  concentrated 
center  load  of  35.5  pounds. 

(2)  Instability  analysis  under  uniform  compression  stress  in  the 
x-direction. 

For  the  deflection  analysis  a  quadrant  of  the  plate  was  idealized 
by  36  elements  as  shown  in  Figure  V-8b.  The  discrete  element  analysis  resulted  in 
a  center  deflection  of  0.0208  inches  corresponding  to  a  center  concentrated  load  of 
35.5  pounds.  From  Reference  18  the  center  deflection  is  determined  by  the  formula 
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(b)  Idealization  of  Plate  Quadrant 


Figure  V-8.  Orthotropic  Rectangular  Plate 
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which  results  in  a  deflection  of  0.025  inch.  It  is  noted  that  the  discrete  element 
approach  gives  a  stiffer  solution.  In  view  of  the  fact  that  the  Hearmon  formula  is  in 
itself  an  approximation,  the  agreement  between  the  two  solutions  is  considered  to  be 
satisfactory. 


The  buckling  stress  of  the  orthotropic,  simply  supported  plate 
was  determined  by  the  expression 


(X 
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From  Reference  19  the  valve  of  c  is  given  as 

Dq 

c  =  (2.0  +  2  -  ...  —  ) 

>x°y 

which  results  in  a  computed  buckling  stress  of  > 


(V-27) 


10,160  lb/in. 


2 


For  the  discrete  element  instability  analysis  the  buckling  mode  is 
assumed  to  be  a  double  buckle,  i.e.,  the  aspect  ratio  2  plate  will  buckle  like  two 
square  plates.  By  considering  antisymmetry  in  the  x-direction  and  symmetry  in  the 
y-direction  it  is  only  necessary  to  analyze  a  quadrant  of  the  plate.  The  idealization 
is  the  same  as  for  the  deflection  analysis  (Figure  V-8b)  and  compressive  edge 
forces  of  1000  Ibs/in.  are  imposed  in  the  x-direction.  Analysis  by  the  subject  pro¬ 
gram  resulted  in  inplane  stresses  of  10,000  lb/in.2  in  the  x-direction  and  a  critical 
inplane  load  factor  (1/1)  of  1.781.  The  predicted  buckling  stress  is  therefore  17,810 
lb/in.2  which  is  considerably  larger  than  the  10,160  lb/in.2  from  Reference  19.  This 
difference  is  consistent  with  the  deflection  analysis  in  that  the  discrete  element 
analysis  again  gave  a  stiffer  result  than  the  alternate  analytical  solution. 
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2.  Analysis  of  Practical  Complex  Conditions 

The  trapezoidal  plate  shown  in  Figure  V-9  furnishes  a  more  practical 
problem  involving  geometric  irregularities.  The  plate  contains  a  centrally  located 
reinforced  hole  and  is  assumed  to  be  simply  supported  on  all  edges.  The  mechanical 
loading  consists  of  a  uniform  lateral  pressure,  inplane  loads  in  both  the  x  and  y 
directions,  and  balancing  shear  flows  along  the  edges.  The  <rx  stresses  are  uni¬ 
form  and  the  o-  stresses  vary  linearly  as  indicated  in  Figure  V-9.  In  addition  the 
plate  is  subjected  to  linear  temperature  gradients  through  the  thickness  which  pro¬ 
duce  thermal  moments. 

Since  the  plate  and  loading  conditions  are  symmetric  about  the  x-axis, 
only  a  half  of  the  plate  is  treated  in  the  analysis.  Figure  V-10  shows  the  half-plate 
idealization  consisting  of  43  elements  and  utilizing  both  quadrilateral  and  triangu¬ 
lar  plates.  Oblique  coordinates  (x',  y ')  were  used  along  the  tapered  edge.  The 
distributed  edge  loads  and  shear  forces  were  transformed  into  equivalent  concen¬ 
trated  x  and  y  forces  at  the  edge  node  points,  and  the  lateral  pressure  was  pro¬ 
rated  as  concentrated  z  forces  to  the  node  points.  The  thermal  moments  were 
hand  computed  from  Equation  V-2. 

The  assumed  material  properties  E  and  a  as  functions  of  temperature 
are  shown  in  the  following  table: 


Temp 
°  F 

E 

106  lb/in.2 

a 

10“6  in./in.-°F 

100 

10.65 

12.6 

200 

10.25 

12.86 

300 

9.82 

13.09 

400 

9.35 

13.32 

500 

8.73 

13.54 

Poisson's  Ratio  was  assumed  constant  at  0.3. 

A  complete  analysis  was  performed  including  inplane,  out-of-plane,  and 
instability  analyses.  The  inplane  direct  stresses  are  shown  in  Figure  V-ll  for 
several  rows  of  elements.  The  stresses  for  the  elements  along  the  x-axis  are 
plotted  in  Figure  V-lla  and  the  stresses  corresponding  to  elements  along  the  sides 
parallel  to  the  y-axis  are  given  in  Figui'es  V-llb  and  V-llc.  The  results  show  that 
the  direct  stresses  are  tension  for  all  the  elements  in  the  reinforcement  around  the 
hole.  Since  the  reinforcement  is  the  cooler  portion  of  the  plate,  the  tension  stresses 
are  reasonable. 

The  out-of-plane  displacements  along  the  ordinates  y  =  0  (x-axis)  and  y=7 
inches  are  plotted  in  Figure  V-12. 
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22  in. 


Section  A-A 

Figure  V-9.  Complex  Trapezoidal  Plate 


Stresses  —  1000  lb/in.  Stresses  —  1000  Ib/in.  Stresses  —  1000  Ib/in. 
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Figure  V-ll.  Inplane  Stresses  for  Trapezoidal  Plate 
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The  inplane  load  factor  (1/  X)  determined  from  the  instability  analysis 

is  1.3  which  gives  critical  stresses  of  c rx  =  1.3(2000)  =  2600  and  cv  =  1.3 

cr  ycr 

(1600)  to  1.3(800)  =  2080  to  1040  or  an  average  <Tycr  °*  1560  Ib/in.  . 


Since  no  test  data  or  alternate  analytical  solution  for  this  plate  is  avail¬ 
able,  it  is  not  possible  to  make  a  comparison  of  the  discrete  element  solutions 
with  other  results.  Approximate  buckling  stresses  can  be  obtained,  however, 
by  assuming  a  20  x  22  inch  simply  supported  rectangular  plate  without  a  hole  and 

loaded  by  uniform  x  and  y  edge  forces.  On  this  basis  the  critical  cr  and  a 
•  x  y 


stresses,  assuming  cr 

2  '2 
lb/in.  and  1450  lb/in. 


0.6  cr  ,  are  determined  from  Reference  20  to  be  2420 
x 

respectively.  These  critical  stresses  are  in  agreement 


with  the  results  from  the  discrete  element  analysis. 
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CHAPTER  VI 

CYLINDER  ANALYSIS  PROGRAM 


A.  SCOPE 

The  cylinder  program  can  be  used  to  analyze  the  structural  behavior,  i.e., 
displacements  and  stresses, for  heated  stiffened  and  unstiffened  fuselage  segments 
and  cylinders.  The  skin  material  may  be  either  orthotropic  or  isotropic  for  an 
untapered  structure  (a  body  of  revolution  with  the  generatrix  parallel  to  the  longi¬ 
tudinal  axis).  For  a  tapered  section,  however,  the  skin  material  must  be  isotropic. 
The  internal  members  (stiffeners  and  frames)  are  permitted  to  be  composed  of  an 
isotropic  material  which  need  not  be  the  same  as  the  skin  material.  The  temperature 
dependence  of  the  material  properties  is  taken  into  account. 

For  illustration  a  simplified  fuselage  section  is  shown  in  Figure  Vl-la.  In 
Figure  VI- lb  is  shown  an  idealization  scheme  which  contains  most  of  the  per¬ 
missible  elements.  For  use  in  idealizing  the  structure,  the  cylinder  program  ac¬ 
commodates  the  following  types  of  elements: 

(1)  Axial  force  elements,  for  idealization  of  longerons  or  longitudinal 
stiffeners. 

(2)  Axial-flexural  element,  employed  in  the  representation  of  frame  elements. 

(3)  Triangular  plate  element,  used  to  represent  either  the  skin  or  the  plate 
portions  of  a  stiffened  bulkhead. 

(4)  Quadrilaterial  plate  element,  for  idealization  of  portions  of  the  skin  and 
bulkheads. 

The  force-displacement  relationships  for  these  elements  were  discussed  in 
Chapter  III. 

The  existing  capabilities  of  the  program  with  respect  to  the  size  of  the 
problem  which  can  be  handled  arc  limited  to  the  inversion  of  a  150th  order  stiffness 
matrix.  This  limitation  cannot  be  precisely  defined  in  terms  of  the  permissable 
number  of  node  points  or  elements;  it  can  only  be  stated  that  no  more  than  150  dis¬ 
placement  degrees  of  freedom  can  remain  after  application  of  the  displacement 
boundary  conditions. 
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Figure  VI- 1.  Fuselage  Segment 
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The  applied  loads  consist  of  concentrated  forces  in  the  x,  y,  z  directions 
and  moments  about  the  x,  y,  z  axes  at  each  reference  point.  If  a  pressure  load  is 
exerted  on  the  shell,  then  it  is  necessary  for  the  analyst  to  redistribute  the  pressure 
load  in  the  form  of  concentrated  loads  at  the  reference  points.  The  program  can 
accommodate  up  to  10  applied  load  conditions  per  structural  problem, 


The  relationships  for  computing  the  axial,  or  inplane,  thermal  forces  are 
coded  in  the  program.  The  computation  of  these  thermal  forces  utilizes  the 
elemental  geometric,  temperature,  and  material  properties  which  are  input  quanti¬ 
ties.  A  single  average  temperature  is  assigned  to  the  axial  and  axial-flexural 
elements.  For  the  skin  elements  (triangular  and  quadrilateral  plates)  the  temper¬ 
atures  are  specified  at  each  reference  point.  The  temperature  of  a  given  skin 
element  is  computed  as  the  simple  average  of  its  corner  point  temperatures. 

The  out-of-plane  thermal  moments  must  be  hand  computed  and  entered  as 
input  to  the  program.  For  the  axial -flexural  element  an  average  thermal  moment 
is  entered  as  part  of  the  element  input  data.  For  the  skin  elements  the  hand 
computed  distributed  thermal  moments  about  the  local  x  and  y  axes  are  stated 
at  each  reference  point.  In  computing  the  concentrated  corner  thermal  moments 
for  a  given  element,  the  distributed  moments  along  the  edges  are  assumed  to  be 
the  average  of  the  distributed  moments  at  the  corners: 


As  in  the  Plate  Analysis  Program,  the  Cylinder  Program  is  capable  of 
accommodating  "oblique'1  support  conditions.  In  some  problems  the  fuselage 
segment  may  be  constrained  to  displace  in  the  direction  of  axes  other  than  those 
employed  in  the  definition  of  the  behavior  of  the  structure  as  a  whole.  Such  con¬ 
ditions  can  be  handled  simply  by  specifying  the  coordinates  of  three  points  which 
define  the  special  axes.  The  program  transforms  the  elemental  relationships 
at  all  the  affected  reference  points  into  equations  referenced  to  the  special  axes, 
and  boundary  conditions  can  be  defined  with  respect  to  the  latter. 


In  addition  to  possessing  the  capabilities  described  above,  the  program  has 
been  designed  to  accommodate  relationships  for  instability  analysis,  two  additional 
types  of  elements,  and  the  capability  to  deal  with  systems  of  approximately  480th 
order.  To  a  limited  extent,  these  additional  capabilities  are  coded  or  contained 
in  the  existing  program  but  were  not  checked  out  as  to  operational  correctness  as 
of  the  conclusion  of  the  subject  study. 
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The  two  above-cited  elements,  which  were  intended  to  complement  elements 
(1) - (4) ,  are  the  following: 

(5)  Composite  plate  element  of  arbitrary  planform,  used  in  idealizing  skin 
and  bulkhead  panels. 

(6)  Composite  axial-flexural  element,  representive  of  longerons,  stiffeners, 
and  frames. 

These  composite  elements,  are  assembled  from  more  basic  elements  in  the 
Plate  Analysis  Program  and  the  One-Dimensional  Element  Analysis  Program,  re¬ 
spectively.  By  incorporating  into  the  present  program  the  appropriate  parts  of 
these  two  programs,  the  elemental  matrices  for  the  composite  elements  are 
computed.  Essentially,  this  process  consists  of  assembling  the  complete  matrices 
for  the  composite  elements  and  reducing  out  the  equations  pertaining  to  the  inter¬ 
mediate  points  which  arc  not  reference  points  on  the  fuselage  structure.  This 
produces  the  element  matrices  referenced  to  the  attachment  points  between  the 
composite  eldment  and  fuselage  structure. 

With  regard  to  planned  capabilities  for  larger  order  systems,  the  program 
has  been  designed  to  accommodate  a  maximum  of  80  reference  points  per  idea¬ 
lization.  ^ince  there  can  be  sLx  degrees  of  freedom  (3  linear  and  3  angular  dis¬ 
placements)  at  each  point,  a  total  of  480  degrees  of  freedom  can  exist  in  an  idea¬ 
lization.  An  important  limitation  is  that  sufficient  "displacement"  and  "force" 
boundary  conditions  must  be  applied  so  that  no  more  than  238  degrees  of  freedom 
remains  in  the  problem.  A  "displacement"  boundary  condition  represents  re¬ 
straint  against  a  given  displacement  component.  A  "force"  boundary  condition 
exists  when  a  force  component  at  a  given  point  is  known  to  have  zero  value  under 
all  load  and  temperature  conditions.  In  cither  case,  the  effect  of  applying  such 
a  condition  is  to  remove  one  equation  and  the  corresponding  unknown  from  the 
problem. 


i 

] 


i 


B.  THEORETICAL  BASIS 

The  theoretical  basis  of  this  program  parallels,  for  the  most  part,  the 
approach  taken  in  the  plate  program  (Chapter  IV).  The  differences  lie  mainly 
in  the  use  of  three,  rather  than  two,  coordinate  dimensions  for  each  element  and 
the  fact  that  in  the  cylinder  analysis  the  membrane  and  bending  behaviors  are 
incorporated  in  a  single  computational  process. 

It  is  of  interest,  however,  to  describe  the  analytical  procedures  involved  in  the 
reduction  of  the  governing  stiffness  equations  by  virtu  re  of  geometric  and  force 
boundary  conditions.  These  procedures  provide  the  theoretical  basis  for  the  planned 
extension  of  the  program  to  accommodate  larger  order  problems  .  In  accordance 
with  the  displacement  approach  to  matrix  structural  analysis,  the  complete  set  of 
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analysis,  the  complete  set  of  force  displacement  equations  are  in  general  given  by: 
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Equation  (VI-1)  has  been  arranged  and  partitioned  as  follows: 

The  first  partition  refers  to  forces  and  displacements  affected  by 
force  boundary  conditions,  tluisjp^j  0. 

The  second  partition  refers  to  forces  and  displacements  unaffected 
by  either  force  boundary  conditions  or  displacement  boundary 
conditions. 

The  third  partition  contains  the  forces  and  displacements  affected 
by  displacement  boundary  conditions  thus  0  and  {^3} 

the  reaction  forces  at  the  support  points. 


(VIrl) 


It  is  possible  Co  reduce  Equation  (VI -1)  so  that  only  the  unaffected  applied 
forces  and  moments  {p,  j  ,  and  displacements  (  1)  appear  in  the  relationship. 

This  is  accomplished  by  simply  removing  the  third  partition  and  operating  upon  the 


remaining  K  and  N  matrices  and  thermal  forces  j  P°  1  as  follows:  The  K 
matrix  is  reduced  by:  ^ 

Huf  hJ  -  K]  hi]"1  h„]  <v>-2> 

tv]  hi]* hi] 
h2)  h,H«„l 

Similar  to  the  reduction  of  the  [k]  matrix,  the  ^  N  j  matrix  is  reduced  by: 

Mrf-  h2]  -  hlhihhz]  m~3> 

The  reduction  of  the  thermal  force  column  is  accomplished  by: 

KJ  BF-h°}  -  [hi!  *  hil]  hh1  K)  <"-» 

The  reduced  form  of  Equation  (VI-1)  can  now  be  written  as 

K}  WllFT  Kl-K}  RF  W-5' 

whore  Wrft  Mhf+Hrf 
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The  displacements  ^A0j,  which  are  unaffected  by  boundary  conditions,  can  be 
solved  for  by  the  expression: “ 


Mrf} 


(VI-6) 


Utilizing  the  known  |  A2  j-displacements,  the  displacements,  -T A ,  associated  with 
the  force  boundary  conditions,  can  be  determined  by  : 


(VI -7) 


Having  determined  all  the  displacements,  the  displacements  of  the  individual 
elements  are  established  and  the  elemental  stresses,  or  stress  resultant  forces, 
are  computed. 


If  critical  stresses  are  to  be  determined  (i.e.,  if  an  elastic  instability  analysis 
is  to  be  performed)  j  P2  }  and  |p°  j  are  set  equal  to  zero,  the  j^Kj  Rp  matrix 

is  reduced,  and  an  equivalent  reduced  ^Nppj  matrix  is  derived.  Theficlpp  matrix 
is  reduced  by  partitioning  it  as  follows:  L  i 
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(VI-8) 


(VI-9) 


An  equivalent  reduced  matrix  is  obtained  by  subtracting  [  Kj^^^from 

the  reduced  matrix  which  is  obtained  in  a  similar  manner  to 

That  is,  first  the  [k  |  matrix  is  partitioned  as  shown: 

RFT 
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(Vl-10) 


and  then  reduced  by  the  expression: 


[k]rftr  [KRFT  J  "  [KRFT  J  [  KRFT  J  [krft  J  (VI  11) 

l.  u  Z  1  11  1 Z 


The  equivalent  reduced  incremental  matrix  is  obtained  by 

Irfr 


Kfr]  [ K  ] rftr  ~  tK  ]rfi 

Lq. 


(VI-12) 


For  the  purpose  of  an  elastic  stability  analysis,  Equation  (VI-1)  may  be 
restated  in  the  reduced  form  as: 


[K]  RFR  T  [NRFlJ 


Eq. 


u  1 
w 


(VI-13) 


Setting 


Fx  I  f  Fxa  1 

>  and  <  pXa  >  equal  to  zero,  the  matrix 

z  j  L <z  J  Ei 


IS 


Eq. 

multiplied  by  the  scalar  X  ,  and  Equation  (VI- 13)  is  rearranged  as  follows: 

-WbfrK 


i- 

X  j  W 


1 


RFRJ 

Eq. 


(VI-14) 


The  scalar  X  ,  which  is  the  eigenvalue,  and  the  eigenvector,  .  (which  is 

the  relative  magnitude  of  the  displacements  and  thus  represents  the  buckled  shape) 
are  determined  through  matrix  iteration. 
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C,  ILLUSTRATIVE  EXAMPLE 

1.  Stress  Analysis  for  Unheated  Conditions 

Figure  VI-2a  illustrates  the  conditions  of  analysis  for  a  ring-stiffened 
cylinder.  The  cylinder  is  cantilevered  and  subjected  to  an  applied  concentrated  load 
at  the  free  end.  The  conditions  shown  were  duplicated  in  a  test  performed  at  the 
NACA  Structures  Laboratory;  results  of  the  test  were  described  in  Reference  21. 

Due  to  symmetry,  only  one  half  the  structure  need  be  considered  in 
analysis.  The  analytical  idealization  sjppears  in  Figure  VI-2a.  Rectangular  plates 
in  plane  stress  (Element  4)  are  employed  in  representation  of  the  skin;  the  flexural 
stiffness  of  these  plates  has  been  neglected.  The  ring  segments  are  axial-flexural 
elements  (Element  2);  their  behavior  is  limited  to  flexure  in  the  plane  of  the  ring, 
shear,  and  direct  axial  loading.  The  middle  line  of  these  elements  coincides  with 
the  middle  surface  of  the  skins  so  that  the  "eccentricity"  is  zero. 

Results  of  the  discrete  element  analysis,  the  solution  obtained  from  beam 
theory,  and  the  test  data  appear  in  Figure  VI--2e.  Only  the  longitudinal  stresses  in 
the  skin  are  compared  but  comparisons  of  other  stress  components  can  be  shown 
to  follow  the  same  trends  and  lead  to  identical  conclusions.  The  beam  theory  results 
are  grossly  in  error  at  ail  points  and,  in  certain  locations,  even  fail  to  predict  the 
correct  sign  (tension  or  compression)  of  the  resulting  stress.  The  accuracy  of  the 
discrete  element  solution,  on  the  other  hand,  is  excellent  at  all  points.  Discrete 
element  solutions  to  this  problem  have  been  published  by  other  investigators  (see 
References  22  and  23).  Generally,  the  latter  have  utilized  shear  panels  in  idealiza¬ 
tion  of  the  skins,  but  their  results  agree  closely  with  the  discrete  element  solution 
of  Figure  VI-2c. 

2.  Cylinder  Thermal  Stress  Analysis 

Anderson  and  Card  (Reference  24)  have  recently  described  elevated  tem¬ 
perature  tests  of  ring-stiffened  cylinders.  One  such  test  specimen  is  shown  in 
Figure  VI-3.  Due  to  the  imposition  of  heat  the  cylinder  skin  assumes  the  longitudinal 
and  circumferential  temperature  profiles  shown  in  Figure  VI-4,  resulting  in  a  state 
of  thermal  stress.  Onty  one-half  the  length  of  the  specimen  is  shown  since  the 
temperatures  are  symmetric  about  ring  number  5.  Also,  there  is  symmetry  about 
the  z-axis  so  that  only  one-quarter  of  the  complete  cylinder  need  be  considered  in 
analysis. 


The  analytical  idealization  is  shown  in  Figure  VI-5.  As  in  the  previous 
example,  the  skins  are  idealized  with  use  of  Element  4  while  Element  2  is  employed 
in  idealization  of  the  ring  segments.  In  the  present  case,  however,  there  is  an 
eccentricity  between  the  middle  line  of  the  rings  and  the  middle  surface  of  the  skin 
and  this  was  taken  into  account  in  the  analysis  performed.  Pertinent  material  pro¬ 
perties,  as  given  in  Reference  24,  are  listed  in  Figure  IV-6. 
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Figure  VI-5.  Idealization  of  Heated  Cylinder 
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The  results  of  analysis  are  shown  in  Figure  VI-7.  This  figure  is  a  re¬ 
production  of  the  one  shown  in  Reference  24  with  the  present  results  being  given  by 
the  heavily  dashed  lines.  As  indicated  the  discrete  element  solution  corresponds 
with  test  data  to  approximately  the  same  extent  as  the  analytical  approach  proposed 
in  Reference  24.  Elementary  theory  is  seen  to  be  entirely  inadequate  for  prediction 
of  the  correct  results. 
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(a)  Variation  of  axial  thermal  stress  in  longitudinal  direction. 


Figure  VI-7.  Comparison  of  Theoretical  and  Experimental 
Thermal  Stresses 
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